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REVISION HISTORY 


This manual is the APMATH64 Manual, Volume 2, 8698-7482-991. The letter 
shown under the revision number column indicates the portion of the 
part number that changes for each revision. The last entry is the 
latest revision to this manual. 


The revision history begins with this manual. 


Deleted Utilities Library, deleted the 
LPSPFI subroutine, added internal subroutine 
information, and added 16 new routines. 


Added new routines to Basic Math Library, 
Doubie Precision Library, and Matrix 
Algebra Accelerated Math Library. 12/87 


NOTE: For revised manuals, a vertical line "|" outside the left 
margin of the text signifies where changes have been made. 


NOTE TO READER 


This is the second volume of the APMATH64 Manual. 
Volume 2 is comprised of part 2 of Appendix A. Note 
that Appendix A continues through Volumes 1, 2, and 
3. The page numbers are listed consecutively through 
the volumes. 


The APMATH64 Manual has three indices located at the 
end of Volume 3 and two at the end of Volume 4. The 
first index (Appendix I) is a list of the APMATH64 
routines in page order by type. The second index 
(Appendix J) is an alphabetical list of all the 
APMATH64 routines. The third index is a key word 
index of the APMATH64 routines. The fourth index 
(Appendix L) is an alphabetical list of the 
APMATH64/MAX routines. The fifth index is a key word 
index of the APMATH64/MAX routines. 
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DESCRIPTION: This routine first calls HTRIDI to reduce A toa 
real symmetric tridiagonal matrix using unitary 

similarity transformations. IMTQL2 is then called to 
determine the eigenvalues and eigenvectors of the 
real tridiagonal matrix. IMTQL2 uses the implicit QL 
method to compute the eigenvaiues and accumulates the 
QL transformations to compute the eigenvectors. 
Finally, HTRIBK is called to backtransform the eigen- 
vectors to those of the original matrix. 


If N is less than or equal to zero, then IERR is set to 
999999. If N is greater than NM, then IERR is set to 
1G*N. If more than 39 iterations are required to 
determine an eigenvalue, the subroutine terminates 

with IERR set equal to the index of the eigenvalue 

for which the failure occurs. In this case, the 
eigenvalues in W should be correct for indices 


ap ie Se ee cne e1genvaiues ase a vermineda wi 


iterations, then IERR is set to zero. 


The function selector, MATZ, may be made functional 
in a future release as follows: If MATZ = 9, then 
only the eigenvalues will be determined; otherwise, 
both the eigenvalues and eigenvectors will be 
determined. 


With the exception of error code 999999 and the 
nonfunctionality of the selector flag, this routine is 
functionally the same as the FORTRAN routine of the 
Same name found in the "Matrix Eigensystem Routines - 
EISPACK Guide", 2nd edition, by B.T. Smith, et al., 
Springer-Verlag (1976). For further information, 

refer to pages 235-239 of the EISPACK Guide. 


The execution time for this routine is highly data 
dependent. 


EXAMPLE: 
Input: 

NM = 4 

N = 4 

AR : 3.8 1.¢ G.2 G.G 
1.0 3.0 G.8 G.G 
g.8 9.0 Le 1.9 
0.0 g.0 ua” | 1.9 
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REARREREREE P REKRREEERER 
* 4 & « 
* RIGRS * -—- REAL SYMMETRIC EIGENSYSTEM SOLVER -—— * EIGRS * 
& & ® & 
REREREREEE REEKEEKERE 
PURPOSE: To determine eigenvalues and eigenvectors of a real 


symmetric matrix 


Ce A =a acto 


CALL FORMAT: CALL EIGRS(NM,N,A,D,E,Z,IERR) 

PARAMETERS : = Integer row dimension of matrices A and Z 

= Integer order of matrix (N .LE. NM) 

Floating-point input matrix 

= Floating-point output vector (eigenvalues) 

= Floating-point scratch vector 

= Floating-point output matrix (eigenvectors) 

IERR = Integer error flag set if routine does not converge 
within 38 iterations (refer to IMTQL2). 


NO YP & 3 
rT) 


NOTE: The dimension of matrices A and Z is NM£*N. 
The dimension of matrices D and E is N. 


DESCRIPTION: EIGRS first reduces the full matrix to tridiagonal 
form by Householder's method, diagonalizing the resulting 
matrix by the QL algorithm (using implicit origin 
shifts). The APAL subroutines used to accomplish this, 
TRED2 and IMTQL2, are based on the FORTRAN programs of 


the same name found in the “Matrix Eigensystem Routines - 
EISPACK Guide" by B.T. Smith et al., Springer-Verlag 
(1976). 
EXAMPLE : 
NM = 5 
N =4 
A bs 5.2 4.9 1.8 1.24 
4.08 5.9 1.9 1.9 
1.80 1.8 4.8 2.0 
1.8 1.89 2.8 4.9 
8.86 8.8 G.f 8.8 


D: 18 2.8 5.8 18.9 
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REREREEREE RREREKEEKE 


& bd = & 
* HTRIBK * —-—- COMPLEX HERMITIAN EIGENVECTORS -—— * HTRIBK * 
x * & 2 
HREREERRER RREREKREEKE 
PURPOSE: To form the eigenvectors of a complex Hermitian matrix, 


A, by back transforming those of the corresponding 
real symmetric tridiagonal matrix determined by the 
routine HTRIDI. 


CALL FORMAT: CALL HTRIBK(NM, N, AR, AI, TAU, M, ZR, 21) 


PARAMETERS : NM = Integer input scalar 

Row dimension of the matrices 

N = Integer input scalar 
Order cf matrix A and column dimen 
the matrices. N must be less than or equal 
to NM. 

AR = Floating-point NM by N input matrix 
The strict lower triangle of the first N rows 
contains information about the unitary trans- 
formations used in the reduction by HTRIDI. 
The remaining elements are ignored. 

AI = Floating-point NM by N input matrix 
The full lower triangle of the first N rows 
contains information about the unitary trans- 
formations used in the reduction by HTRIDI. 
The remaining elements are ignored. 

TAU = Floating-point 2 by N input matrix 
Contains the remaining information about the 
unitary transformations. 

M = Integer input scalar 
Number of eigenvectors to be back transformed. 

ZR = Floating-point NM by N input/output matrix 
On input, the columns of ZR contain the eigen- 
vectors to be back transformed in their first N 
elements. On output, the first M columns and N 
rows contain the real parts of the transformed 
eigenvectors. 

ZI = Floating-point NM by N output matrix 
The first M columns and N rows contain the 
imaginary parts of the transformed 
eigenvectors. 
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RRERRECKRAE RREREKAEEE 

& & * * 

* HTRIDI * -——— COMPLEX HERMITIAN TRIDIAGONALIZATION ——- * HTRIDI * 

* & 7 x 

RERRARRTERE REEKREEAEEE 

PURPOSE: To reduce a complex Hermitian matrix, A, to a real 
symmetric tridiagonal matrix using unitary similarity 
transformations. 


CALL FORMAT: CALL HTRIDI(NM, N, AR, AI, D, E, E2, TAU) 


PARAMETERS : NM = Integer input scalar 
Row dimension of the matrices 

N = Integer input scalar 
Order of matrix A and column dimension of 
the matrices . N must be less than or equal 
to NM. 

AR = Floating-point NM by N input/output matrix 
On input, the first N rows of AR contain the 
real parts of the elements of A. The last 
NM - N rows are ignored. Only the full lower 
triangle of AR need be supplied. On output, 
the strict lower triangle of AR contains 
information about the unitary transformations 
used in the reduction. The full upper 
triangle of AR is unaltered. 

AI = Floating-point NM by N input/output matrix 
On input, the first N rows of AI contain the 
imaginary parts of the elements of A. The 
last NM - N rows are ignored. Only the strict 
lower triangle of AI need be supplied. On 
output, the full lower triangle of AI contains 
information about the unitary transformations 
used in the reduction. The strict upper 
triangle of AI is unaltered. 

D = Floating-point output vector of length N 
Contains the diagonal elements of the 
tridiagonal matrix. 

E = Floating-point output vector of length N 
Contains the subdiagonal elements of the 
tridiagonal matrix in its last N-1 elements. 
The element E(1) is set to zero. 

E2 = Floating-point output vector of length N 
Contains the squares of the corresponding 
elements of E. 

TAU = Floating-point 2 by N output matrix 
Contains the remaining information about the 


unitary transformations. 
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| REREERAEEE RARRERRERER 
|* * * * 
| * IMTQL1 * -——— DIAGONALIZE TRIDIAGONAL MATRIX --- * IMTOL1 * 
|* * * * 
j**teeeeeeee RRAKKEREKKE 
| 

| 

| PURPOSE: To determine the eigenvalues of an N by N 

real symmetric tridiagonal matrix using the 


implicit QL method. 


CALL FORMAT: CALL IMTQLI (N, D, E, IERR) 


PARAMETERS : N Integer input order of the matrix 

D = Floating-point input/output vector 
Vector of length N containing the diagonal 
elements of the symmetric matrix on input; 
vector of length N containing the eigenvalues 
on output. 

E = Floating-point input vector 
Vector of length N containing the subdiagonal 
elements of the symmetric matrix. The 
subdiagonal is contained in elements E(2) 
through E(N); E(1) is arbitrary. 

IERR = Integer output error status 


IERR = 9: No errors encountered, normal 
completion. 
IERR = -l: The routine received an invalid 


| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| input argument, N < l. 

| TERR > @: The routine was unable to finish 
| because more than 39 iterations 
| were required to determine an 

| eigenvalue. IERR is set to the 
| index of the offending eigenvalue. 
| The eigenvalues in D are correct 
| for all preceding indices, but are 
| unordered. 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 

| 


DESCRIPTION: IMTQL1 determines the eigenvalues of a symmetric 
tridiagonal matrix using the QL algorithm with 
implicit origin shifts at each iteration. 


Upon convergence, the eigenvalues are ordered in 
ascending order. 


The vector E is destroyed by this routine. 


IMTQL1 is based on the FORTRAN program found in 
the EISPACK GUIDE, 2nd ed.; B.T. Smith, et al.,; 
Springer-Verlag, 1976. That program in turn is 
based on an Algol procedure discussed by Martin 
and Wilkinson, NUM. MATH.,12, 1968, pg. 377. 
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* IMTQL2 * 


2 & 
Fe SESELE SE & 1 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


To determine e 
symmetric tridiagonal matri 
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RRRERERERE 
2 = 


--- DIAGONALIZE A TRIDIAGONAL MATRIX ---— * IMTQL2 * 


# + 
kkkeekeeee 


CALL IMTQL2(NM,N,D,E,Z,IERR) 


NM = Integer row dimension of matrices A and 2 
N = Integer order of matrix (N .LE. NM) 
D = Floating-point input/output vector 


Diagonal elements on input; 

Eigenvalues in ascending order on output 
Floating-point input vector 

Codiagonal elements 

Z = Floating-point input/output matrix 

For eigenvectors of sym.tridiag. matrix: 
Nth-order identity matrix on input; 
Eigenvectors on output 

For eigenvectors of full sys. matrix: 
Trans.matrix from TRED2 on input; 
Eigenvectors on output 
Integer index of eigenvalue if convergence not 


bi 
MW 


IERR = 
obtained by 38 iterations, else g 

NOTE: The dimension of arrays A and Z is NM*N, 
The dimension of arrays D and E is N. 


IMTQL2 diagonalizes an N-by-N tridiagonal matrix 

using the implicit QL algorithm (Martin and Wilkinson, 
Num. Math. 12, 377(1968); Dubrulle, Num. Math. 15, 459g 
(1978)). The initial diagonal begins at D(1l), and the 
codiagonal at E(2). At each iteration, a new tridiagonal 
matrix is formed, is stored by overwriting the previous 
result, and continues until convergence, or 39 iterations 
have passed. If convergence does not occur by 34 
iterations, IERR is set equal to the index of the sought 
eigenvalue, and the routine is exited. Previously 
calculated results are valid. The transformation 
matrices are accumulated and the results stored in column 
Order in matrix Z. 
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REREERAAEE RAERRRKEEEE 
& ® *® * 
* RS * -— REAL SYMMETRIC EIGENSYSTEM SOLVER ——- * RS * 
2 2 * x 
RVERATERER RRARRREERET 
PURPOSE: To determine the eigenvalues and eigenvector of a 


real symmetric matrix, A. 


CALL FORMAT: CALL RS(NM, N, A, W, MATZ, Z, FV1, FV2, IERR) 


PARAMETERS : NM = Integer input scalar 
Number of rows of matrices A and Z 
N = Integer input scalar 


Order of matrix A and column dimension of | 
matrices A and Z. N must be less than or equal 
to NM. 

A = Floating-point NM by N input matrix 
The first N rows contain the matrix and the 
last NM ~ N rows are ignored. Only the full 


lower triangle of the matrix need be supplied. 


W = Floating-point output vector of length N 
Contains the eigenvalues of A in ascending 
order. 


MATZ = Integer input scalar 
MATZ is not currently used. 
Z = Floating-point NM by N output matrix 
The first N elements of the j-th column of 2 
is the eigenvector that corresponds to the 
j-th eigenvalue in W. The last NM-N 
elements in each column are not altered. 
FV1l = Floating-point work area vector of length N 
FV2 = Floating-point work area vector of length N 
IERR = Integer output scalar 
Error code as described below. 


DESCRIPTION: This routine first calls TRED2 to reduce A toa 
symmetric tridiagonal matrix using and accumulating 
orthogonal similarity transformations. IMTQL2 is 
then called to determine: the eigenvalues and eigen- 
vectors of the original matrix from the symmetric 
tridiagonal matrix. IMTQL2 uses the implicit QL 
method to compute the eigenvalues and accumulates the 
QL transformations to compute the eigenvectors. 


If N is less than or equal to zero, then IERR is set to 
999999. If N is greater than NM, then IERR is set to 
19*N. If more than 38 iterations are required to 
determine an eigenvalue, the subroutine terminates with 
IERR set equal to the index of the eigenvalue for which 
the failure occurs. 
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RRERREATERE REREREEKRER 
& & & x 
* SIMPLE * ~———- REVISED SIMPLEX ——- * SIMPLE * 
. ® * x 
RRAERRRERRE REREREKEEE 
PURPOSE: To solve a linear programming problem that is in the 


standard form: 
maximize Z = C'* X 


subject to A*xX=B 


and X(j) >= 8, for j = 1 to N 
where B(i) >= 8, for i=1toM 
CALL FORMAT: CALL SIMPLE(M,N,MP2,NP1,KI,NS,S,IRN,ICP,B,C,WRK, 
X,Y,Z,IB,KO) 
PARAMETERS : M = Integer input scalar 
Number of constraints (rows in A). 
N = Integer input scalar 


Number of variables (columns in A). 
MP2 = Integer input scalar 


MP2 = M + 2 

NP1 = Integer input scalar 
NPl=N+t+1 

KI = Integer input vector of length 19 
Contains the program control parameters. If 


any of these parameters is less than or equal 
to zero, then a default value is supplied for 
that parameter. The parameters are: 
KI(1) = Input basis flag. KI(1) > @ indicates 
that an initial basis is supplied in 
IB. Default = No initial basis. 
KI(2) = Iteration limit. Default = 4 * N+ 19g 
KI(3) = Inversion interval. Default = M/2 + 5 
KI(4) = Zero tolerance exponent. The zero 
tolerance value = 9.5 ** KI(4). 
Default = 2g. 
KI(5) = Partial pricing step size. 
Default = min (N, max(29,N/29)). 
NOTE: The default value is also used 
if KI(5) > N and a value of 29 is 
used if @ < KI(5) < 29. 
KI(6) to KI(1%) are reserved for future use. 


NS = Integer input scalar 
Number of nonzero elements in A. 

Ss = Floating-point input array of length NS 
Contains the nonzero elements of A stored by 
columns. 


IRN = Integer input array of length NS 
Contains the row numbers (in A) that 
correspond to the nonzero elements in S. 
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The problem must be stated in the standard form: 
maximize Z = C'* X 


subject to A* xX =B8B 
and X(j) >= 9, for j = 1toN 


where B(i) >= 9, for i=1toM 


Therefore, it is the res 


Therefore ; onsibility of the user to: 

(a) Convert a minimization problem to a maximization 
problem by replacing C with -C. 

(b) Convert inequality constraints to equality con- 
straints by adding a slack variable or sub- 
tracting a surplus variable. 

(c) Ensure that B(i) >= @ by multiply the i-th con- 
straint by -1.9 if B(i) < @. 

(dj) Ensure that the decision variables are con- 
strained to be nonnegative. If X(j) is uncon- 
strained in sign then replace it by the 
difference of two new nonnegative variables. 


In this variation of the two phase, revised simplex 
method, a composite problem is formed (virtually) in 
SIMPLE that includes both the actual (phase 2) 
objective equation and the artificial (phase 1) 
objective equation as constraints making a total of 
MP2 constraints. The variables for the internal com- 
posite problem are: 

X(1) to X(N) - The actual decision variables 
X(N+1) - The artificial objective 

X(N+2) to X(N+M+1) - The artificial variables 
where X(N+l+i) is the artificial variable for the 
i-th constraint. 


The variables X(%@) and X(N+tl) to X(N+M+1) are virtual 
variables and, thus, do not use any storage space. 


X(@) must always be a basic variable and IB(l) must 
always be zero. X(N+1l) must be a basic variable 
during phase one and IB(2) must equal N+l whenever 
X(N+1) is basic. At least one artificial variable 
(including X(N+1l)) must always be basic. During 
phase two, any artificial variables in the basis will 
have a value of zero. Generally, during phase two, 
only one artificial variable will be basic and it 
will be X(N+1); however, this need not be the case. 
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EXAMPLE: 


Given a problem in standard form where 


AS dia tac Bec de, Bie Bs 
G9. O G@. 3. 41. 2. 
Ze Sa “Ba: 25. Wes, Be 
G9. -@s: Bz: Be. -Se Be 
33 BO. De BO. Gs 3. 
the inputs are: 
M = 5 
N = 13 
MP2 = 7 
NP1 = 14 
KI Py g, g, Q, g 
NS = 22 
S. “@°2s9-2sy. Bey espe) 
Sag aue Zar 347 Lev 
IRN : 1, 3, 5S- Il, 3, 
Sr ay 4, 4, 5, 
ICP : 1, 4, 6, 8, 19, 
22, 23 
B 2 1l4.5 252% 21s% 38; 
Cc : 9 eae se bas 
Bir Bey Ber Ge 
IB : Don't care since KI(1l) = 9g 
The outputs are: 
xX Py g., 53:7 g., E Per g., 
B., @. 
Y 
Z = 177.9 
IB by g, 14, 9, 2s 8, 
KO bY g, 7+ Ty Gg, Ty 


G. @. 
@. 3B. 
2. @. 
2. 3. 
g. l. 
= ee ae 
Ley Lig 
l, 4, 
l, 2, 
12, 14, 
34. 

7. 6 


8., O., 


l. 
g. 
B. 
g. 
g. 


1g., 
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@. @. G. 
1. @. @. 
g. 1. g. 
Oe Be he 
G9. GO. . 
2s oie Sy 
Levy das 

3, 2, 4, 
4, 5 

18, 19, 29, 
we. Baus a9 
4., g., Ge, 


21, 


3 2.0000, 9.6667, 3.9999, 1.4815, 1.5556, y6, y7 
where y6 and y7 not of interest (scratch) 
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EXAMPLE: 


4.988 1.909 9.82808 7.989 
@.988 -3.9009 8.908 -2.9008 -7.998 


TY 
ht 
uo 
RM 
i} 
® 


A( INPUT) 


A( OUTPUT) 13.998 3.923 8.877 9.889 4.796 
2.898 -8.765 1.786 -9.425 -%.778 


V(INPUT) = 8.289 1.899 2.999 2.5282 2.22828 
MAXA = a 2 4 5 8 ll 
NN = 5 
MA = 3 

- NWA = 19 
KKK = 3 
V(OUTPUT) = -%.943 9.563 G@.245 9.483 8.315 
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e 
° 


| EXAMPLE 


WARRA R tas] las] tax] 
SS e- " a 
WMWMeBwmW ws as} ta] R 
i Ba N a oO 
BMeBR®RA RW ta] ta] 
NOM WN WR wn Py =f 
SS Ss tas) ax] tax} 
MRNA OM Be et =f Ve) 
t i q qo 
WM BMWA TR ~® Lax] 
HAaN at & N tas] tos) 
in un 1 ] 
now 
N 
z & i aA fx] ica 
oe 5 
cu) sa 
3 uy 
Q, 43 
G =) 
H Oo 


- 257 


Page A 


FPS 868-7482-991C 


APPENDIX A 


RAREEEREREE RERREREEEE 
. 7 & * & 
* VASORT * -——— VECTOR SORT ALGEBRAIC VALUES ——— * VASORT * 
& . x = 
RERRERAAEE REEKREREELE 
PURPOSE: To sort a vector into an ascending vector 


of algebraic values using Quicksort. 
CALL FORMAT: CALL VASORT(A,I,N,W) 
PARAMETERS : = Floating-point vector to be sorted in place 
= Integer element step for A 
= Integer element count 


= Floating-point vector of at most 2*log2(N) words 
of contiguous space for working stack of 


i - 


i Oe we Te ae 


DESCRIPTION: VASORT sorts elements of a vector into an ascending 
vector of algebraic values by the method of 
Quicksort (Hoare's partition-exchange sort) in 
place. The procedure iteratively partitions 
the vector creating two subvectors, one whose 
values are less than or equai to the value 
initially at the middle location, and the other with 
elements greater than or equal to that value. 

This chosen value ends up in its true 
(post-sorted) position between the two subvectors. 
The half-way location was chosen for initial trial 
comparison in order to speed the sort when the 
Original vector is already partly ordered. 


After each partition, first and last locations of 
the larger subvector are stored in a pointer stack, 
which can accumulate no more than log2(N) pairs, 
and the process of partitioning is continued on the 
smaller subvector. The process of comparison and 
partitioning is continued until no subvectors 
remain. The vector is then completely sorted. 


EXAMPLE: 
N=5 
A(input) : 5.8 4.8 -9.5 -1.9 8.9 
A(output) : -1.9 -0.5 4.8 5.8 8.9 
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* VSORT * 
# # 
ERALALRARE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


APPENDIX A 


SERRARRERRER 
x 2 
-——- VECTOR SORT WITH INDICES -—- * VSORT * 
& & 
REREREREKE 


To sort a vector into an ascending vector of 
algebraic values using Quicksort. When the elements 
of the A vector are swapped, corresponding elements 
of the P vector are also swapped. Typical use of 

the P vector is to record the original indices of the 
sorted vector. 


CALL VSORT(A,I,P,J,N) 


= Floating-point vector to be sorted in place 
= Integer element step for A 

Integer or real vector of starting indices 
Integer element step for P 

Integer element count 


24 GH YP 
iu 


VSORT sorts elements of a vector into an ascending 
vector of algebraic values by the method of 
Quicksort (Hoare's partition-exchange sort) in 
place. The procedure iteratively partitions 

the vector creating two subvectors, one whose 
values are less than or equal to the value 
initially at the middle location, and the other with 
elements greater than or equal to that value. 

This chosen value ends up in its true 
(post-sorted) position between the two subvectors. 
The half-way location was chosen for initial trial 
comparison in order to speed the sort when the 
Original vector is already partly ordered. 


After each partition, first and last locations of 
the larger subvector are stored in a pointer stack, 
which can accumulate no more than log2(N) pairs, 
and the process of partitioning is continued on the 
smaller subvector. The process of comparison and 
partitioning is continued until no subvectors 
remain. The vector is then completely sorted. 
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RRERRAEETREE RERRERERELN 
. . od & 
* ACORT * -——- AUTO-CORRELATION (TIME-DOMAIN) --—— * ACORT * 
x . 7 & x 
RREARARREKER REREEKEXAE 
PURPOSE: To perform an auto-correlation operation on a vector 


using time-domain techniques. 
CALL FORMAT: CALL ACORT(A,C,N,M) 


PARAMETERS : A = Floating-point input vector 

= Floating-point output vector 

N = Integer element count for C 
(Number of lags) 

M = Integer element count for A 

p 


io) 
i 


Pat a+ Fo ok ee 1 nt 


(nN ~ ama “a ~~ mene ol mmm Ss 
(NV VEVelLUL SLCMWSlLLS YeeUpPy GCULSSCLULSL 


addresses.) 


DESCRIPTION: C(m)=SUM(A(m+q-1)*A(q)), 
for gq=l1 to M-m+l 
m=l1 to N 


ACORT uses time-domain techniques (compare with ACORF) 
to compute the auto-correlation function. This routine 
needs less storage than ACORF, and runs faster when N 
and/or M is small. The resultant vector C must not 
overlay the source vector A. 


EXAMPLE : 
N = 3 
M=5 
A: 1.9 2.8 3.8 4.0 5.8 
C : 55.8 489.89 26.9 
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RARREEEEEEE 
® ® 


* BLKMAN * 
2 ® 
KREKEAERERE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


EXAMPLE: 


——— BLACKMAN WINDOW MULTIPLY -——— 


To multiply a vector by a Blackman window. 
CALL BLKMAN(A,1I,C,K,N) 


= Floating-point input vector 

Integer element step for A 

= Floating-point output vector 

= Integer element step for C 

= Integer element count (a power of 2) 


ZRAAH YP 


rT /ATAN 


ffmiTie/ de I/N)) 


g Z*Cco VCmmay ye” 
+8 .98 


for m=l1 to N 


Multiplies the elements of the vector A by 
and 
N must be 


an N element Blackman window function, 
stores the results in the vector C. 
a power of 2. 


I=l 

K=1 

N = 8 

A: 1.8 1.8 1.2 1.0 
1.9 1.9 1.2 1.9 

C: 8.998 8.866 9.340 G.774 
1.9808 8.774 8.348 9.866 
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THLINC = Floating-point input scalar containing 


the phase increment threshold (used to 
obtain more confident phase estimates near 
sharp zeros) 


THLCON = Floating-point input scalar containing 


the phase consistency threshoid 

= Integer work area vector of length 39 
used for various software stacks during phase 
unwrapping 


IXCXST = Integer input scalar X and CX input status: 


@ if X is provided as input and. 
CX is not provided as input 
1 if X is not provided as input and 
CX is provided as input 
2 if both X and CX are provided as input 


IAUXST = Integer input scalar AUX input status: 


@ if AUX is not provided as input 
1 if AUX is provided as input 


th 


TIPHWST = Integer inout scalar chase unwr in nle 


ree ad et Gel hee Co do pst unwrapping only 
status: 
@ if complex cepstrum is desired 


1 if phase unwrapping only is desired 


NOTE: For APFTN64 calls to CCEPS, the dimension 


DESCRIPTION: 


1) 


2) 


7) 
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of arrays X, CX, and AUX must be greater 
than or equai to NFFT2 and the dimension 
of array WMD must be greater than or equal 
39. 


See "Programs for Digital Signal Processing", 


IEEE Press, 1979. 


Input parameters are checked for out of range 
conditions. If any errors are detected, then 
SSUC gets the appropriate error code (2.9 - 9.9) 
and CCEPS returns. 

If IXCXST=9 then X is used to compute CX. 

If IXCXST=1 then CX is used to compute X. 

Note that in this case the vector X will occupy 
NFFT2 words in Main Memory but only the first NX 
elements of X will be used in further calculations. 
If IAUXST=9 then X is used to compute AUX. 

Each of the NFFT2 elements of CX and AUX are 
divided by 2.9 to match IEEE formulation. 

If the first element of CX is less than 9.9 

then SNX = -1.9 

else SNX = +1.9 . 

The magnitude of the spectrum is computed and 
stored in the real positions of AUX; 

the phase derivative of the spectrum is computed 
and stored in the imaginary positions of AUX; 
and twice the linear phase estimate (mean of the 
phase derivative) is computed for use in the 
phase unwrapping computation. 


CX(OUT) : —-1.6639 

8.9543 

3.5771 

AUX (OUT) : @.8359 

6.7434 

1279 .4929 
SNX = 1.8889 
SPX = -2.1909 
ssuc = 9.9908 
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%.G88S 
1.4985 
6.8888 


-2.6149 
-2.7728 


~2.9325 


-5.9134 
3.9149 


6.8889 


415.6262 
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* CCORT * 
* * 


RRERBRRRRRRE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


EXAMPLE: 


APPENDIX A 


RERERELEEE 
* * 
——— CROSS-CORRELATION (TIME-DOMAIN) ——— * CCORT +# 
* & 
RERKRRARRE 


To perform a cross-correlation operation on two 
vectors using time-domain techniques. 


CALL CCORT(A,B,C,N,M) 


= Floating-point input vector (operand) 

= Floating-point input vector (operator) 
Floating-point output vector 

= Integer element count for C (number of lags) 
= Integer element count for A and B 

(Note vector elements occupy consecutive 


addresses.) 


220A YP 
u 


C(m)=SUM(A(m+q-1)*B(q)); 
for q=l1 to M-m+l1 
and m=1 to N 


CCORT uses time-domain techniques (compare with CCORF) to 
compute the cross-correlation function. This routine 
needs less storage than CCORF, and runs faster when N 
and/or M is small. 


N = 3 

M= 4 

A: 1.2 2.0 3.2 4.9 
B: 169.9 29.0 308.0 48.9 
C : 398.8 2009.8 118.9 


FPS 8698-7482-981C Page A — 273 


APPENDIX A 


RRERERERRE REREKRERERE 


bd * & * 
* COHER * -——— COHERENCE FUNCTION ——— * COHER * 
& . x . 
ERREKEAEAVEE BPTTERKIETIEE 
PURPOSE: To compute the coherence function, given the 


auto-spectra of two signals and the cross-spectrum 
between them. 


CALL FORMAT: CALL COHER(A,B,C,D,N) 


PARAMETERS : A = Floating-point input vector 

(Auto-spectrum) 

B = Floating-point input vector 
(Auto-spectrum) 

C = Complex-floating-point input vector 
(Cross~spectrum) 

D = Floating-point output vector 
(Coherence function) 

N = Integer element count 
(Note vector elements occupy consecutive 
addresses.) 


DESCRIPTION: D(m)= (R(C(m))**2+I(C(m))**2)/(A(m)*B(m)); for m=l1 to N 


EXAMPLE: 


a 
ul 
Ww 


90ND Y 


se 08 ee ee 
_ 
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RRRAEAEREEETE RREEEEKEREE 
x bd & 5 4 
* DECFIR * —-- DECIMATION ——— * DECFIR * 
. 7 & & ed 
RRERRERERE REERREKRERE 
PURPOSE: To FIR filter an input vector using a convolution 


technique incorporating decimation by a factor D. 
Typically, the input vector is a digital signal 
requiring low pass filtering and the operator vector 
is the array of pre-determined filter coefficients. 


CALL FORMAT: CALL DECFIR(A,B,C,D,N,M) 
PARAMETERS : 


= Floating-point input (undecimated) vector 
= Floating-point input operator vector 


BPilanati ng-point Outout vector 


eww Beg wep vueweve 


= Integer input decimation factor (D > @) 
= Integer input element count expected for C 
when convolving without decimation 
(NOTE: the actual size of the output vector 
C will be [(N-1)/D]+1) 
M = Integer input element count for B 
(NOTE: element count for A must be N+M-1) 


20OQO OW YP 
i} 


DESCRIPTION: C(m) = SUM (A(D*(m-1)+q) * B(q) for q=l to M 
and m=l1 to [(N-1)/D]+l 
(NOTE: This assumes that the operator array B 
is loaded with the elements arranged in 
reverse order. Thus: 
B(1) = Mth operator point 
B(2) = (M-1L)th operator point 


e 


B(M) lst operator point) 

For references see: 

(1) A. Peled and B. Liu, “Digital Signal Processing: 
Theory, Design and Implementation." John Wiley, 
1976. 

(2) R. E. Crochiere and L. R. Rabiner, "Optimum FIR 
digital filter implementation for decimation, 
interpolation, and narrow band filtering," IEEE 
Trans. Acoust. Speech Signal Processing, 
vol ASSP~23 pp 444-456, Oct. 1975. 


This routine performs a convolution on the decimated 


operand A with the operator B. The results are 
stored in {(N-1)/D]+l elements of vector C. 
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RRREKEUTERE RAEREEREEEE 
& * * & 
* ENVEL * -—— ENVELOPE DETECTOR -—— * ENVEL * 
« . 2 & 
RERARERLREX RREKEREKER 
PURPOSE: To obtain the envelope of a vector X(t). 


CALL FORMAT: CALL ENVEL(X,E,N) 


PARAMETERS : Xx 
E 
N 


Floating-point input vector 
Floating-point output envelope vector 
Integer element count (a power of 2) 


DESCRIPTION: E(t) = SORT { x(t)**2 + H{X(t)}**2 } for t=l to N 
where: H{xX(t)} = Hilbert transform of X(t). 


For references see any standard text on 
communication theory, viz. “Communications Systems 
and Techniques," M.Schwartz, W.Bennet, & Stein, 
McGraw Hill. 

This routine starts by obtaining the Hilbert 


transform of the input vector. The formula shown 
above is then used to extract the envelope. 


EXAMPLE: 


Xs 2.7 1.6 8.3 4.2 9.7 14.1 3.6 G.5 


E: 2.72 4.32 8.82 4.38 11.33 14.73 9.21 9.85 


FPS 868-7482-991C Page A - 279 
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* HANN * 


RRETLATEZE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


EXAMPLE: 


APPENDIX A 


RREREREAREE 


——- HANNING WINDOW MULTIPLY -———- * 


To multiply a vector by a Hanning window. 
CALL HANN(A,I,C,K,N,F) 


= Floating-point input vector 

Integer element step for A 

= Floating-point output vector 

= Integer element step for C 

= Integer element count (a power of 2) 

= Integer normalization flag 
@ for unnormalized Hanning window 
(peak window value=1.9) 
1 for normalized Hanning window 
(peak window value=1.63) 


MARA p 


N should be a power of 2. I£ not, HANN sets N to the 
next lower power of 2. For further information see 
Digital Time Series Analyses, Otnes and Enochsen, Jonn 
Wiley '72, page 294. 


C(m)=W*A(m) *(1.9-COS(2*PI*(m-1)/N); for m=l1 to N 


where: 

W g.5 for F=9 

W = 9.8165 for F=l1 

N= 4 

F= 9 

A:1.8 1.8 1.89 1.9 

C: 8.8 8.5 1.8 @.5 

N= 4 

Fe=l 

A: 1.98 1.98 1.989 1.989 
C : G.08 g.82 1.63 9.82 
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RERREARREX 
* ® 


* HLBRT * 
* * 
HeRREREARE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


EXAMPLE: 


APPENDIX A 


RAERERERAETE 
* * 
——— HILBERT TRANSFORMER ——- * HLBRT * 
& & 
REREREEAEE 


To obtain the Hilbert transform of an analytic 
signal. 


CALL HLBRT(X,H,N) 
Floating-point input vector 


Floating-point output Hilbert transformed vector 
Integer element count (a power of 2) 


am x 
wou 


F{a#{x(t) = -J * F{x(t)} for t=1 to N 

where: F{X(t)} = Fourier transform of X(t). 
H{X(t)} = Hilbert transform of X(t). 
di = SQRT(-1) 


(1) A real to complex FFT of X(t) is obtained. 

(2) Real components of the result are multiplied by 
el 

(3) Positions of the real and imaginary components 
are switched. 

(4) A complex to real inverse FFT is performed on the 
results of step 3. 
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EXAMPLE : 
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de 


= 
MPmRwuoan naw fF WN & 


— 


ty 
w 
u 


= 128 


ig 


: 18 * SIN(i * 19 * 2*PI/128) + 
28 * SIN(i * 28 * 2*PI/128) + 


306 * SIN(i * 38 * 2*PI/128) 


for i 
RC(3j) 


-9.2847 
9.8183 
-8.52908 
8.5493 


%.18262 


FobhWVVuay 


-9.2451 
-8.8955 
@.3127 
8.4627 
8.3854 


1 to 128 


A(j) 


1.9909 
-9.8897 
1.9494 
-$.25989 


~ 129A 


erases 


9.2885 
@.1145 
9.8661 
9.1981 
9.1478 
9.39854 


AL(j) 


89599 .9 
82335.4 
27198.6 
19844.6 


TAGS. 4 


dee “ER WS te Oo 


13563.2 
12748.7 
12632.4 
11397 6L 
8957.3 
8121.7 


R(j) 


89599.9 
25512.8 
-69112.8 
~37858.5 


299%2 a 


theme te Ss 


24552.1 
-38917.1 
-18999 .6 

35262.2 

1979759 
=52577%9 
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EXAMPLE: 


A: 8.8 18.0 20.0 3.09 4.8 58.8 6.89 78.9 


B: 26.0 3.8 58.0 6.9 
C: 3.86 4.80 6.8 7.9 
R = 4.8 
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Inverse transform: 
x(k) = suM { X((r-1)*df-Pl*df) * 
EXP(j*2*pi*(r-1)*d£*(k-1)*dt) } 
for r = 1 to NF 
where dt = 1/XM. 


Thus the same formula used for the forward transform 
may be used for the inverse transform if here 

W = -2*pi*(k-1)/XM and Fl and NF replace Tl and 

NT respectively. If the r=1 component X(1) is 
input, it must have an imaginary part equal to @. 


The DFT is produced by the modified Goertzel algorithm 

as described in 

(1) A.V. Oppenheim and Schafer, “Digital Signal 
Processing," Prentice Hall, 1975 

and 

(2) F. Bonzanigo, "An improvement of Tr 
phase uUnvrapping algorithm,” IEEE T 
Feb. 1978, pp. 194-195 

Additionally, an exponential factor has been used 

to account for any offset of the input values from 

zero (Tl or Fl). 


Inverse times are approximately double for forward 
times after the NT and NF values are interchanged. 


EXAMPLE: 

Fl = 9.9 

Tl =1.9 

NT = 8 

NF = 4 

XM = 8.9 

I =1 

A(INPUT) : 1.0 9.8 -1.0 9.8 1.8 G9.8@ -1.0 Gf 


B(OUTPUT) : (8.9, 8.9)(9.8, 9.8)(9.9,-4.9) (9.8, 8.8) 


B(INPUT) : (4.9, &.8)(0.8, @.@) 


A(OUTPUT) : 8.8 9.8 -8.09 9.8 8.0 G.8 -8.8 9.8 
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& z= 


* RFTIIT * 
* * 
Seeeeeeees 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


EXAMPLE: 


APPENDIX A 


REREERREEE 
* = 


-——- REAL FFT WITH QUARTER INTERPOLATION -——- * RFTII * 


® *& 
SREREKERESR 


To perform an in-place real-to-complex forward or 

a complex-to-real inverse fast Fourier transform (FFT) 
including the case of N=64K via quarter interpolation 
in the 4K cosine table. 


CALL RFTII(C,N,F) 


C = Floating-point input/output vector 
N = Integer input element count (power of 2) 
F = Integer input direction flag: 

+1 for forward 

-l1 for inverse 


See RFFT. 

N= 4 

F = 1 (Forward) 

C(IN) 3: 1g.8 12.9 18.9 12.8 
C(OUT) : (89.9,9.f) ( 8.9,89.f8) 

N= 4 


F = ~l (Inverse) 


: (89.09,09.f8) ( G.9,09.f8) 
C(OUT) : 88.0 89.08 88.8 89.0 
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REREEEKEEE RERRERAKEE 
7 eo z 2 
* TCONV * -—-— POST-—TAPERED CONVOLUTION (CORRELATION) ——- * TCONV * 
2 ® : x & 
RERZEREKEE BERKREREREE 
PURPOSE: To perform a post-tapered convolution or correlation 


operation on two vectors. 


CALL FORMAT: CALL TCONV(A,I,B,J,C,K,N,M,L) for correlation 
CALL TCONV(A,1I,B(N),J,C,K,N,M,L) for convolution 

PARAMETERS : = Floating-point input vector (operand) 

= Integer element step for A (>) 

= Floating-point input vector (operator) 

= Integer element step for B (<% => Convolution) 


= Integer element step for C 

= Integer element count for C 
= Integer element count for B 
= Integer element count for A 


CE WZRAGWH D 
id 
2 
> 
] 
+ 
3 
°] 
] 
5 
3 
+ 
2 
= 
cy 
{= 
ct 
< 
D 
7 
- 
3 
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FORMULA: C(m)=SUM(A(m+tq-1)*B(q)); 
for q=1 to R 
and m=l1 to N 


where: 
R=MIN(M, L-M+1) 


DESCRIPTION: TCONV performs either a correlation (I and J positive) or 
a convolution (I positive and J negative) operation 
between the L-element operand (trace) vector A and the 
M-element operator (kernel) vector B. The N-element 
result vector is stored in C. TCONV automatically 
inserts zeros into the calculation if N+M-1 exceeds the 
operand length L, thus saving storage and zeroing of 
N+M-1-L extra operand elements. (Compare with CONV.) 


EXAMPLE: 


C4 
" 
N 


3.8 5. 


w 


11.9 19.9 
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RUERREREEETE RRREKKRERE 
* J x & 
* TRANS * ——— TRANSFER FUNCTION —— * TRANS * 
> x rd & 
REKREREREE REREKRERERE 
PURPOSE: To perform a complex transfer function calculation by 
dividing the cross-spectrum by the auto-spectrum, 


CALL FORMAT: CALL TRANS(A,B,C,N) 


PARAMETERS : A = Floating-point input vector 

(Auto~spectrum) 

B = Complex-floating-point input vector 
(Cross-spectrum) 

C = Complex-floating-point output vector 
(Transfer function) 

N = Integer element count 
(Note vector elements occupy consecutive 
addresses.) 


DESCRIPTION: R(C(m))+I(C(m))=(R(B(m))+I(B(m)))/A(m); for m=l to N 


EXAMPLE: 
N = 3 
A: 1.2 2.8 3.8 
B: (1.9,2.8) (3.9,4.8) (5.8,6.9) 
C : (1.8,2.9) (1.5,2.8) (1.67,2.9) 
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RUERAEAEEREE RERARKRERETE 
z *® 2 x 
* VAVLIN * ——-— VECTOR LINEAR AVERAGING --—— * VAVLIN * 
® 7 & J 
HRBVRRRRRRKE REKTARERERE 
PURPOSE: To update the linear average of a sequence of vectors 


to include a new vector. 


CALL FORMAT: CALL VAVLIN(A,I,B,C,K,N) 


PARAMETERS : A Floating-point input vector 
I = Integer element step for A 
B = Floating-point input scalar 
(Number of vectors included in current average) 
= Floating-point input/output vector 


Integer alamant step for nr 


amdils ae Woe de aut 44 


= Integer element count 


2m 
i 


DESCRIPTION: C(m)=C(m)*B/(B+1.3) + A(m)/(B+1.9); for m=l1 to N 


EXAMPLE: 
N= 5 
A : 5.8909 18.008 29.980 25.999 39.290 
B : 5.809 
C( INPUT) 2:19.90 18.900 19.980 198.900 18 .GG8¢ 


C(OUTPUT) : 9.167 18.909 11.667 12.599 13.333 
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RERRAETAEEE RERRERERAE 
2 * * & 
* VxCS * -——-— VECTOR MULTIPLIED BY SIN AND COS -— * vxcS #* 
* * (TABLE LOOKUP) * * 
RERRRAARER VERRREEEEN 
PURPOSE: To multiply a vector with the sine and cosine 


s 
of a linearly increasing argument with a given 


ee a esse i Seve ats = 


initial phase. 
CALL FORMAT: CALL VXCS(A,C,K,F,P,N) 


PARAMETERS: A = Floating-point input vector to be multiplied by 
the sine and cosine functions 
C = Complex floating-point output vector 
K = Integer input element step for C 
(K >= 2) 
F = Floating-point input scalar frequency 
P = Floating-point input scalar phase at t= 
= Floating-point output scalar initial phase value 
for next frame 
N = Integer element count 


DESCRIPTION: Re(C(m)) 
Im(C(m) ) 


A(m) * COS((m-1)*F+P) 
A(m) * SIN((m-1)*F+P) 
for m= 1toN 


NOTE: The argumen 


This routine multiplies vector A with a sine and 
cosine function defined by frequency F and initial 
phase P. Straight ROM table lookup is used for 
generating the sine and cosine values and thus this 
routine has limited precision. The initial phase 
value for the next frame is returned in P. 

NOTE: K should be greater than or equal to 2 

sO as not to destroy part of the resultant vector 

C as it is generated. 


EXAMPLE : 


a 'UWAR 
I fT) 
Ow Qn 
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RRERRERERE REKRREREKEE 
* & ed x 
* WIENER * ——- WIENER LEVINSON ALGORITHM -—— * WIENER * 
x & & x 
SERERESAEER SRAELKUITKALECSE 
PURPOSE: To solve a system of single channel normal 


equations which arise in least squares filtering 
and prediction problems. 


CALL FORMAT: CALL WIENER(LR,R,G,F,A,ISW,IERR) 


PARAMETERS : LR = Integer filter length 

R = Floating-point input vector (Auto-correlation 
coefficients) 

G = Floating-point input vector (Cross 
correlation) 

F = Floating-point output vector (Filter 
weighting coefficients) 

A = Floating-point output vector 


(Prediction error operator) 
ISW = Integer input (algorithm switch) 
@ = spike deconvolution 
1 = general deconvoiution 
IERR = Integer output scalar (failure flag) 


DESCRIPTION: WIENER solves: 
1. The following set of LR equations for F; 


SUM [F(p)*R(m-pt1)=G(m); 
for p=l1 to LR and m=1 to LR 


2. The following set of LR equations for A; 


SUM [A(p)*R(m-p+l)=V*D; 
for p=l1 to LR and m=1 to LR 


where, A(1)=1.9 
D=1.9 when m=1 
D=9.8 when m not = 1 
V=A(L)*R(L)+...+A( LR) *R( LR) 
R(-i)=R(i) 


If the algorithm is successful IERR is set to @; 


else it is set to the pass number at which the 
failure occurred. 
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IMAGE PROCESSING LIBRARY 
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a = 
1 


Aan 
Am 


-l (Inverse) 


F 


C( IN) 


( 4.8,8.8) ( 0.9,8.8) ( 8.0,8.8) ( 0.8,08.8) 


(16.8,0.8) (16.8,8.8) (16.8,0.8) (16.9,8.8) 


C(OUT) 


- 395 


Page A 
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IR = Integer input scalar flag: 
“non-zero for correlation 
9 for convolution 


DESCRIPTION: C((itIC-1),(j+JC-1))= 
SUM(A( (itIA+k-2-irbias), (j+JA+1-2-icbias) )*B(k,1) ) 

where i=l to M 

j=l to N 
for k=l to MB 

l=1 to NB 

and IBL=MB*NB-IB1+1 for convolution 
icbias=(IB1-1)/MB 
irbias=(IB1-1)-MB*icbias 
(row and column biases are from the 
initial B(1,1) position. 


CONV2D correlates or convolves a two-dimensional 
operand submatrix A' of A with a two-dimensional 
operator matrix B, and stores the result in 
submatrix C' of C. A one-to-one correspondence 
exists between the elements of A' and C'. 


This routine does not do boundary testing. 
Therefore care must be taken when choosing values 
for IA, JA, and IBl for given values of M, N, MB, 
NB, and IR to avoid using data outside of A when 
computing C'. 


EXAMPLE: 

MA = 9 

IA =1 

JA =1 

M = 7 

N 7 

MB = 3 

NB = 3 

MC = 9 

IC 1 

Jc =l 

A: 6.09 6.89 8.6 6.09 G68 O.8 G8 86.8 G8 
9.9 1.09 1.80 1.0 4.60 4.09 8.9 8.8 8.8 
6.68 1.6 1.86 1.0 4.0 4.09 8.9 86.09 G8.8 
9.8 1.6 1.0 1.6 4.0 4.0 8.089 86.8 G.8 
6.68 1.60 1.6 1.86 4.09 4.0 8.89 6.9 G.8 
8.8 8.89 9.89 6.80 8.89 G09 86.09 OH B.E 
6.68 8.86 89.86 0.9 9.9 G89 8.09 B.GF GG 
9.2 8.2 §.8 €.82 2.8 2.8 2.9 8.2 8.2 
8.8 G8 8.80 6.09 O88 OCH G8 GC.H GB. 
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Here the operator, B, is positioned for processing the initial 
point in a’. 
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Here the operator, 8B, is positioned for processing the initial 
point in A’. 
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Figure A-Z Convoiution 
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EXAMPLE 


ONAN ON N OO /O 
“on tf f ft tw toon 


et U 
Sash ee 


9.09 G09 G9 G8 G8 GG GH G.G 


A 


1801.6 1.86 189 1.8 6.8 


1.9 


G.8 


18 18 8.0 


1.2 


1.8 1.8 
9.6 1.86 1.6 2.86 2.6 1.89 1.89 8.8 


9.6 1.86 1.6 2.8 2.0 1.86 1.89 O.8 


G.8 1.9 
9.9 1.9 


1.8 9.0 


1. 


1.8 1.9 


1.g 


9.9 O96 8.89 G86 6.89 G9 G8 BH 


U 
U 
U 
U 
U 
U 


2.9 
338 
3.2 


3.9 
1.9 


3.0 


3.8 3.9 
1.869 2.8 2.6 


2.8 
3.8 


U 
U 
U 
U 
U 
U 


3.8 2.86 2.8 2.89 2.9 
3.86 2.6 2.9 2.8 2.9 


3.6 1.80 2.8 2.9 


as) 
3.8 


3 


1.9 


2.8 3.86 3.0 3.8 3.80 2.9 


(U indicates unchanged elements of C) 


= Sil 
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This routine differs from GRAD2D in that it can 
perform testing for image boundaries, substituting 
zeros for values that are needed outside the 
boundary. The routine runs somewhat more slowly 
than GRAD2D. 


If testing is employed, zeros are substituted for 
those elements in the formula which fall outside 

of A. This is useful in preventing wrap-around 

and incorrect processing of the columns and rows on 
the borders of A. However, the testing adds pro- 
cesSing time and is unnecessary when there is a 
border of width one around A' which lies totally 
within A. 


If boundary testing is not employed (i.e. B = J) and 
if a boundary of A' coincides with all or part of a 
boundary of A, then boundary effects will be observed 


* : 
mn thea f2FAmnutatian af Ct Tn + = 
in the computation of C'. In the cases of JA=1 


JA+N-1=NA these boundary effects may not be pre- 
dictable since data stored adjacent to A may not be 
predictable. 


ar 
wt oe 


EXAMPLE : 

MA = 8 

NA 8 

IA al 

JA 1 

Cc 64 

MC = 8 

NC 8 

Ic =l 

JC 1 

M = 8 

N = 8 

B =l1 

A: 1.89 1.89 1.89 1.8 1.9 1.9 1.89 1.9 
18 #18 1.2 L:8 1s 1.8 1.8 1.9 
1841.8 1.86 1.8 1821.8 «18 1.9 
1.8 1.89 1.68 2.8 2.8 1.9 1.7 1.9 
1.8 1.8 1.8 2.89 2.8 1.89 1.2 1.8 
1891.8 1.86 1.8 #1686 «168 1.8 1268 
1861.8 1.8 18 1.86 $18 1.8 1.8 
bee. cig Tae tees “Te Lee. 168 Led 
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RARERAREAEEE REREREEERE 
x 2 * 2 
* LAPL2D * ——— LAPLACIAN FILTER -— * LAPL2D * 
& & & x. 
RAERKREAREKEE ReRERRKREEE 
PURPOSE: To filter images for edge enhancement by 

applying a two-dimensional Laplacian 


operator. 
CALL FORMAT: CALL LAPL2D(A,MA,NA,IA,JA,C,MC,NC,IC,JC,M,N,IX) 


PARAMETERS : A = Floating-point input matrix 

(column ordered) 

MA = Integer number of rows of A 

NA = Integer number of columns of A 

IA = Integer initial row of the submatrix A‘* of A 
to be processed (1 < or = IA < or = MA) 

JA = Integer initial column of the submatrix A‘' of A 
to be processed (1 < or = JA < or = NA) 

C = Floating-point output matrix 
(column ordered) 

MC = Integer number of rows of C 

NC = Integer number of columns of C 

IC = Integer initial row of C which locates the. 
submatrix C', where C' will be the processed A' 
(1 < or = IC < or = MC) 

JC = Integer initial column of C which locates the 


submatrix C' (1 < or = JC < or = NC) 
M = Integer number of rows in A‘ 

(1 < or = M < or = MA) 
N = Integer number of columns in A' 


(lL < or = N < or = NA) 
IX = Integer distance to filter side from center of 
Square: side S=2*(IX+1); filter area = S**2 


DESCRIPTION: C'(p,q)= 128 -4*A'(p,q)+A' (p-IX,q)+A' (pt+IX,q) 
+A'(p,q-IX)+A' (p,qtIX) 


Each of the elements in C' is calculated according 
to the above formula, which adds to a bias of 128 
a weighted combination of each pixel and its 4 
horizontal and vertical neighbors at distance IX. 


If a boundary of A' coincides with all or part of a 
boundary of A, then boundary effects will be observed 
in the computation of C'. In the cases of JA<=IX or 
JA+N-IX>=NA these boundary effects may not be pre- 
dictable since data stored adjacent to A may not be 
predictable. Boundary effects will be predictable 

if A' is initially ringed with a known constant, 

such as zero. 
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RERRERAEAREE REEREAERERE 
. . 2 , 
* LPL2DB * -—— LAPLACIAN FILTER WITH BOUNDARY TEST —-—- * LPL2DB * 
* 2 . & 
RREKKEKKKAE BERRETREEE 
PURPOSE: To filter images for edge enhancement by applying 


a two-dimensional Laplacian operator. This 
routine does special boundary testing. 


CALL FORMAT: CALL LPL2DB(A,MA,NA,IA,JA,C,MC,NC,IC,JC,M,N,IX,B) 


PARAMETERS ; A = Floating-point input matrix 
(column ordered) 
MA = Integer number of rows of A 
NA = Integer number of columns of A 


IA = Integer initial row of the submatrix A' of A 
to be processed (1 < or = IA < or = MA) 


JA = Integer initial column of the submatrix A' of A 
to be processed (1 < or = JA < or = NA) 

C = Floating-point output matrix 
(column ordered) 

MC = Integer number of rows of C 

NC = Integer number of columns of C 

IC = Integer initial row of C which locates the 
submatrix C', where C' will be the processed A' 
(1 < or = IC < or = MC) 

JC = Integer initial column of C which locates the 
submatrix C' (1 < or = JC < or = NC) 


M = Integer number of rows in A' 
(1 < or = M < or = MA) 
N = Integer number of columns in A' 


(1 < or = N < or = NA) 

IX = Integer distance to filter side from center of 
Square: side S=2*(IX+1); filter area = S**2 

B = Integer input scalar which is 8 if no boundary 
testing is desired; if not = @, values needed 
outside of A are evaluated as zeros 


DESCRIPTION: C'(p,q)= 128 -4*A'(p,q)+A' (p-IX,q)+A' (p+IX,q) 
+A'(p,q-IX)+A'(p,qtIX) 


Each of the elements in C‘ is calculated according 
to the above formula, which adds to a bias of 128 
a weighted combination of each pixel and its 4 
horizontal and vertical neighbors at distance IX. 


This routine differs from LAPL2D in that it can 
perform testing for image boundaries, substituting 
zeros for values that are needed outside the 
boundary. The routine runs somewhat more slowly 
than LAPL2D. 
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REKKRERERE . REEREREKTEE 
* * & z 
* MED2D * --- MEDIAN FILTER -— * MED2D * 
x * * x 
RERAREARAEAE REEAREREETEE 
PURPOSE: To filter out noise in images by replacing 


each pixel with the median value of the pixels 


— soe we NG we eRe ¥ NG ate Ww biwe 


in a square window centered around the pixel. 
CALL FORMAT: CALL MED2D(A,MA,IA,JA,C,MC,IC,JC,M,N,IX,H,L) 


PARAMETERS: A = Floating-point input matrix 
(column ordered) 
MA = Integer number of rows of A 
NA = Number of columns of A) 
A = Integer initiai row of the submatrix A‘* of A 
to be processed (1 < or = IA < or = MA) 
JA = Integer initial column of the submatrix A' of A 
to be processed (1 < or = JA < or = NA) 
C = Floating-point output matrix 
(column ordered) 
MC = Integer number of rows of C 
(NC = Number of columns of C) 
IC = Integer initial zow of C which locates the 
submatrix C', where C' will be the processed A' 
(1 < or = IC < or = MC) 
JC = Integer initial column of C which locates the 


submatrix C' (1 < or = JO < or = NC) 
M = Integer number of rows in A' 

(1 < or = M < or = MA) 
N = Integer number of columns in A' 


(1 < or = N < or = NA) 

IX = Integer distance to median filter side from 
center of square: side S=(2*IX)+1; 
filter area = S**2; IxX>g 

H = Floating-point vector histogram used as a work 
area 

L = Integer input scalar length of H = 
2**(number of bits per pixel) 


DESCRIPTION: C'(p,q)=median of all elements A'(t,u), 
p-IxX<=t<=ptIx, q-IX<=u<=q+Ix 


For each of the elements in A' a histogram is 
formed from the median of the elements within 
+ or ~ IX row and column distance from the element. 


The median is found via a fast algorithm pu 
in: 
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REREREEERE RERREAKREEE 


® * x 7 
* MOVREP * ~——— SUB-IMAGE MOVE AND LEVEL REPLACE --- * MOVREP * 
* * * x 
PES SESE E SE | REEKEREREE 
PURPOSE: To simply move a sub-image A' of an image A and/or 


to replace each pixel value with another value 

as specified in the lookup table, vector T, whose 
elements are the new values and whose subscripts 
are the original pixel values + l. 


CALL FORMAT: CALL MOVREP(A,MA,IA,JA,C,MC,IC,JC,M,N,T,NT) 


PARAMETERS : A Floating-point input matrix 
(column ordered) 


= 
MA = Integer number of rows of A 


(NA = Number of columns of A) 


IA = Integer initial row of the submatrix A' of A 
to be processed (1 < or= IA < or = MA) 

JA = Integer initial column of the submatrix A' of A 
to be processed (1 < or= JA < or = NA) 

C = Floating-point output matrix 
(column ordered) 

MC = Integer number of rows of C 


(NC = Number of columns of C) 

IC = Integer initial row of C which locates the 
submatrix C’' of C; C' will be the processed A' 
(1 < or= IC < or = MC) 


JC = Integer initial column of C which locates the 
submatrix C’ of C (1 < or= JC < or = NC) 

M = Integer number of rows in A’ 
(1 < or = M < or = MA) 

N = Integer number of columns in A' 
(1 < or = N < or = NA) 

T = Floating-point input vector pixel replacement 
table 

NT = Integer input scalar length of vector T = 


2**(# Of bits per pixel) 
(NT = @ indicates only submatrix move is 
desired) 


DESCRIPTION: For pixel replacement, 
C*'(p,q)=T(FIX(A' (p,q) )+1) 


For submatrix move, 
C'(p,q)=A' (p,q) 
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* RFFT2D * 
* * 
REREKKEREE 


PURPOSE: 
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RERRARKEEEE 

‘ * *& 

-——- REAL TO COMPLEX 2—DIMENSIONAL FFT —— * RFFT2D * 
= & 

RERRREREERE 


To perform an in-place two-dimensional real-to- 
complex forward or a complex-to-real inverse fast 


Fourier transform (FFT). 


CALL FORMAT: CALL RFFT2D(C,N1,N2,F) 


PARAMETERS : 


C = Floating-point input/output matrix 
(column ordered) 
Nl = Integer number of rows = 
number of real elements per column 
(power of Z2 < or = 16384) 
N2 = Integer number of columns = 
number of real elements per row 
(power of 2 < or = 16384) 
NOTE: N1*N2 must be < or = available main data 
F = Integer direction flag: 
+1 for forward 
-l1 for inverse 


DESCRIPTION: Forward: RFFT2D performs a two-dimensional real to 


complex forward FFT on the Nl by N2 real array C, 
storing the (N1/2 + 1) by (N2/2 + 1) complex array 
result in a special packed complex array form 


ae bd rw +n tt te Zz te 


occupying the same Nl by N2 locations of array C: 


Let El = N1/2 and E2 = N2/2 


R(1,1) 
R(El+1,1) 
R(2,1) 
I(2,1) 


R(E1,1) 
I(El,1) 


R(1,E2+1) R(1,2) I(1,2) ... R(1,E2) I(1,E2) 
R(El1+1,E2+1) R(E1+1,2) I(E1+1,2).. R(E1+1,E2)1I(E1+1,E2) 
R(2,2) R(2,3) R(2,4) ... R(2,N2-1) R(2,N2) 
I(2,2) I(2,3) R(2,4) ... I(2,N2-1) I(2,N2) 
R(E1,2) R(E1,3) R(E1,4) ... R(E1L,N2-1)I(E1,N2) 
I(El1,2) I(E1,3) I(E1,4) ... I(El,N2-1)I(E1,N2) 


The results of a two-dimensional real-to-complex 
forward FFT should be multiplied by 1/(2*N1*N2) for 
proper scaling. 
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LINPACK BLAS LIBRARY 
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REEREKEREER REREREEEEE 
4 *# & & 
* CAXPYN * ——- NESTED COMPLEX A * X + ¥ -— * CAXPYN * 
2 & & x 
REREEEAREEE RERZREREEEE 
PURPOSE: To add a scalar multiple of one complex floating-point 


vector to another complex floating-point vector N times, 
each time for a different pair of vectors and a different 
complex floating-point scalar. The first vector is a 
subset of the vector X, and the second is a subset of the 


vector Y. The scalar is an element of the vector A. 
CALL FORMAT: CALL CAXPYN(ISW,N,M,A,IAO,X,IXI,IX0O,Y,IYI,IYO) 


PARAMETERS : ISW = Integer input scalar. ISW is a function 
selector switch and is treated as a bit 
string with the bits numbered from the 
least significant bit (bit @). If a given 
bit is set (equal to 1), then the function 
option that corresponds to that bit is selected. 
All options are independent of each other and 
are summarized below. 

Bit J: Negate A * X. 

Bit 1: Not used. 

Bit 2: Use conjugate of A. 

Bit 3: Use conjugate of X. 
All other bits are ignored. 


N = Integer input scalar. Number of A* K+ Y¥ 
Operations, i.e., outer loop count. 

M = Integer input scalar. Number of elements in 
each A * X + Y operation, i.e., inner loop 
count. 

A = Complex floating point input vector. Array of 
scalars. 

IAO = Integer input scalar. Outer loop element 
increment for A. 

x = Complex floating point input vector. First 
input vector. 

IXI = Integer input scalar. Inner loop element 
increment for X. 

IxXO = Integer input scalar. Outer loop element 
increment for X. 

Y = Complex floating point input/output vector. 
Second input vector on input. Output vector on 
output. 

IYI = Integer input scalar. Inner loop element 
increment for Y. 

IYO = Integer input scalar. Outer loop element 


increment for Y. 
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Output: Y 
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Hew ow 
Wr Re FR WN RM 


( 3.9,-1.9) 


( 2.9, 9.8) 


( 2.9, 1.9) 


(-1.9, 9.9) 


( 2.8, 8.8) 
( G.8,-2.8) 


(-1.9, 1.2) 
(-2.9,-2.9) 
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SKEREEREEKEE RERAEREREEE 
* *& * * 
* cpoTC * ——— COMPLEX INNER PRODUCT -—— * CDOTC +* 
x & * rd 
RREERRREEEE Pee EE LE SS 
PURPOSE: To sum conjugates of first complex vector 

times elements of second complex vector. 


CALL FORMAT: CW = CDOTC(N,CX,I,CY,J) 


PARAMETERS : N = Integer element count 
CX = Pirst complex floating-point input vector 
I = Integer element step for CX 
CY = Second complex floating-point input vector 
J = Integer element step for CY 

Ch Compiex fioating-point output vaiue 


C 
= 
) 


DESCRIPTION: CW = SUM((R(CX(m))-I(CX(m)))*(R(CY(m))+I(CY(m)))); 
for m=l to N 


CW = (9.9,8.8) if N<l. 


EXAMPLE: 
N = 2 
I=l 
J=l 
CX (9.38,0.48) (0.99,1.92) 


CY : (8.38,8.48) (8.98,9.92) 
: (9.25,-8.9) 
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Z = Complex floating point input/output vector. 
An input only if bit 1 of ISW is set. 
IZO = Integer input scalar. Element increment 


for Z. 
DESCRIPTION: 2Z(jz) = r * Z(jz) + s * SUM[ X(ix) * Y(iy), i=1,M] j=1,N 


wheres ix = (j-l) * IXO + (i-1l) * IXI +1 


iy = (j-1) * IYO + (i-1) * IYI +1 

32 = (4=—1) * TT7N + 7 

je ‘joe sala ~ 

s = 1.8, if ISw({d] = @ 
= -1.9, if ISw(g] = 1 

r= 9.8, if ISsw{l] = g 
= 1.9, if Isw{l} = 1 

X = X , if ISsw(2] = g 
= Conjg(X), if ISw{Z] = i 

¥ = Y , if Isw(3] = @ 
= Conjg(Y), if ISwW(3] = 1 

Z = Zz, if Isw([4] = g 


= Conjg(Z), if ISW[4}] = 1 


and ISW(k]}] = bit k of ISW. 


NOTES: If IZ0 = 8, then CDOTN will set 2Z(1) equal to 
the accumulated sum of all N dot products. If 
ISW{1L] = 1 alsc, then input 2Z({1)}) will be added 


to this sum. 

Memory words occupied by X may intersect those 
occupied by Y. In fact, X and Y may coincide. 
However, memory occupied by Z should not, in 
general, intersect that occupied by X or Y. 


If N < 1, CDOTN returns with no action taken. 


If M < 1 and ISW{1l] = 1, CDOTN returns with 
no action taken. 


If M < 1 and ISW{1l] = 9, CDOTN returns with 
Z(j) = @.@ for j = 1 toN. 


In general, M < 1 implies a zero sum of 
products. 
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* cpoTU * 
x * 
RRERAEEEEEE 


PURPOSE: To 


CALL FORMAT: CW 


PARAMETERS : N 


DESCRIPTION: CW 


EXAMPLE: 


a 


CX 
CY. 
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RERKKEKEEE 
* * 
—— COMPLEX DOT PRODUCT --- * CDOTU * 
*& * 
RHRRERERKE 


compute the inner (unconjugated) product 
two complex vectors. 


= CDOTU(N,CX,I,CY,J) 


= Integer element count 

= First complex floating-point input vector 
= Integer step for CX 

= Second complex floating-point input vector 
= Integer step for CY 


= Complex floating-point scalar output result 
= SUM(CX(m)*CY(m)); for m=1 to N 


= (9.9,9.8) if N<l. 


te 


(8.38,09.49) (8.99,1.99) 
(9.38,-.48) (8.99,9.99) 
(-8.75,8.92) 
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KeEKREKERER RREEKEREREL 
7 ® & 
* CSCAL * -—— COMPLEX SCALING --— * CSCAL * 
2 x * & 
RRRREEKEEE PeeS EL ESSE YS 
PURPOSE: To multiply each component of a vector 


by a complex scalar. 
CALL FORMAT: CALL CSCAL(N,CA,CX,TI) 
PARAMETERS : N = Integer element count 


CA = Complex floating-point scalar multiple 
CX = Complex floating-point input/output vector 


I = Integer step increment for CX 
DESCRIPTION: CX(m) = CA*CX{(m); for m=1 to N 


EXAMPLE: 
N = 3 
I=l 
CA : ( G.2, 1.2) 
CX(INPUT) : ( 1.8, 2.8) ( 3.8, 4.8) ( 5.8, 6.9) 
CX(OUTPUT) : (-2.9, 1.8) (-4.9, 3.8) (-6.9, 5.@) 
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* CSSCAL * 


x * 
REREREAEEE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


EXAMPLE: 
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RERAKEEREE 
& * 


* CSSCAL * 
£ & 


-——— REAL TIMES COMPLEXES -——-— 


RAEEAKKRRRREK 


To multiply the elements of a complex vector 
haev =a raal anal ar 
wy a Car VwVaLah o 


CALL CSSCAL(N,SA,CX,1) 


N = Integer element count for CX 

SA = Floating-point input scalar multiple 

CX = Complex floating-point input/output vector 
I = Integer element step increment for CX 
CX(m) = SA*CX(m); for m=1 to N 

N = 3 

I=l 

SA : 8.5 

CX(INPUT) : (2.9,4.9) (6.9,8.8) (9.9,1.@) 
CX(OUTPUT) : (1.9,2.8) (3.8,4.8) (9.8,98.5) 
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* ICAMAX * 
& & 
SERAELAKAAER 


PURPOSE : 


CALL FORMAT: 


PARAMETERS : 


EXAMPLE: 
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* * 


* ICAMAX * 
® x 


-——- INDEX OF LARGEST COMPLEX ELEMENT ——— 


RRAERRKEREEE 


To calculate the index of the complex 
element of largest real plus imaginary magnitude. 


IMAX = ICAMAX(N,CX,TI) 

N = Integer element count 

CX = Complex floating-point input vector 

I = Integer step increment for CX 

IMAX = Integer value of index with largest components 


cmag(CX(IMAX)) = 
where cmag(C) = 


MAX(cmaq(CX(m)); m=l1 for N 
ABS(R(C))+ABS(I(C)), 


with 1 < = IMAX < = WN. I£ N < 1, IMAX = @. 
N = 3 

I=l 

CX : ( 3.9, 3.8) ( 5.9,-9.9) ( 9.9,13.2) 
IMAX : 2 
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* SASUM * 
* * 
RHERRRRARRE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


EXAMPLE: 
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RRAREKEEEE 

& & 

——- SUM OF MAGNITUDES ——- * SASUM * 
. ® & 
ERRAEREREX 


To sum magnitudes of elements of a real 
vector. 


SW = SASUM(N,SX,1I) 


N = Integer element count 

SX = Floating-point source vector 

I = Integer incremental step for SX 
SW = Floating-point scalar result 


SW = SUM(ABS(SX(m))}); for m=l to N 
N = 3 

SX :-1.89 9.8 5.9 

SW: 6.9 
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REKERREREE RETRAIN 
7 . 2 x 
* SAXPYN * ——- NESTED REAL A * X + Y ———- * SAXPYN * 
x = e 2 & 
RERREERERE RREEREREKEX 
PURPOSE: To add a scalar multiple of one floating-point vector 


to another floating-point vector N times, each time 

for a different pair of vectors and a different scalar. 
The first vector is a subset of the vector X, and the 
second vector is a subset of the vector Y. The scalar 
is an element of the vector A. 


CALL FORMAT: CALL SAXPYN(ISW,N,M,A,IAO,X,IXI,IxX0,Y,1YI,IYO) 


PARAMETERS: ISW = Integer input scalar. ISW is a function 
selector switch and is treated as a bit 
string with the bits numbered from the 
least significant bit (bit 8). If a given 
bit is set (equal to 1), then the function 
option that corresponds to that bit is selected. 


Only bit @ is used in SAXPYN. 


Bit 9: Negate the product term A * X 
before adding to Y. That is, 
compute - A * X + Y instead of 
A * X + Y¥. 

All other bits are ignored. 
N = Integer input scalar. Number of A* X + ¥Y 


operations, i.e., outer loop count. 


M = Integer input scalar. Number of elements in 
each A * X + Y operation, i.e., inner loop 
count. 

A = Floating point input vector. Array of 
scalars. 

IAO = Integer input scalar. Outer loop element 
increment for A. 

Xx = Floating point input vector. First input 
vector. 

IXI = Integer input scalar. Inner loop element 
increment for X. 

IXO = Integer input scalar. Outer loop element 
increment for xX. 

6 = Floating point input/output vector. Second 
input vector on input. Output vector on 
output. 

IYI = Integer input scalar. Inner loop element 
increment for Y. 

IYO = Integer input scalar. Outer loop element 


increment for Y. 
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A : 3.8 -1.9 2.9 

xX =: 24 3.0 

Y +: 7.8 6.9 2.0 3.8 5.2 6.9 
Output: Y :13.8 15.9 G.8 9.9 9.9 12.¢ 
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x x 
* SCNRM2 * 
x = 
SRAREEELITEST 


PURPOSE: To 


CALL FORMAT: SW 


PARAMETERS : N 


SW 


we] 
(= 
VP) 
©) 
og 
tH 
tg 
| 
tI 
o) 
a 
Fs 


EXAMPLE: 


CX 
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REERREKEEX 
x x 
——— COMPLEX EUCLIDEAN NORM —- * SCNRM2 * 
& 7 
PRRBAEVREEEETE 


compute the square root of sum of squares 
elements of a complex floating-point vector. 


= SCNRM2(N,CX,1) 


= Integer element count 

= Complex floating-point input vector 
= Integer step increment 

= Floating-point scalar output result 


| oie 


(9.8,3.9) (4.8,8.f) 
5.0 
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RRREREEREA REKREREREX 
& * x 7 
* spor * —— DOT PRODUCT OF REAL VECTORS -— * SDOT * 
& . *& & 
RRAKEAETRETE RRERKKEERAE 
PURPOSE: To compute the inner (dot) product 


of two vectors. 
CALL FORMAT: SW = SDOT(N,SX,1I,SY,J) 
PARAMETERS: N = Integer element count for SX and SY 
SX = Floating-point input vector 
I = Integer element step for SX 
SY = Floating-point input vector 
J = Integer element step for SY 
SW = Floating-point output value 


eow Ss ieedind! ~dinadind aA a3 


DESCRIPTION: SW=SUM(SX(m)*SY(m)); for m=1 to N 


SW=9.08 if N<l. 


EXAMPLE: 
N = 3 
SX : 1.8 2.8 3.9 
SY : 4.8 2.5 9.8 
SW: 5.2 
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DESCRIPTION: Z(jz) = cr * Z(jz) + s * SUM({ X(ix) * Y(iy), i=1,M] j=1,N 


wheres ix (j-1) * IXO + (i-1) * IXI +1 
SL) YO id) ee 


jz = (j-l) * I1Z0 +1 


| oad 
a) 
' 


s = 1.8, if Isw(d] = g 
= -1.9, if ISW(@9] = 1 
r = §.8, if ISW(1] = @ 
= 1.9, if ISW{lJ] = 1 
and ISW[k] = bit k of ISW. 
NOTES: If IZ0 = @, then SDOTN will set 2Z(1) equal to 


the accumulated sum of all N dot products. If 
ISW(1] = 1 also, then input Z(1) will be added 
to this sum. 


Memory words occupied by X may intersect those 
occupied by Y. In fact, X and Y may coincide. 
However, memory occupied by Z should not, in 
general, intersect that occupied by X or Y. For 
sample applications, see Sections D.4.9 and D.4.11. 


If N < 1, SDOTN returns with no action taken. 


If M < 1 and ISW[{1]} = 1, SDOTN returns with no 
action taken. 


If M < 1 and ISW{1] = 9, SDOTN returns with 
Z(j) = @.@ for j = 1toN. 


In general, M < 1 implies a zero sum of products. 


EXAMPLE: 


Input: ISW = 


g 
N = 2 
M = 3 
IXI 2 
IxO = 1 
Iyt 1 
TYO = @ 
IzO = 1 


X : 3.9 2.8 -1.8 18 8.8 -2.8 


Output: Z : 1.8 -2.@ 
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RERREAEEETE RRRERERERTE 
x ® 7 & 
* SROT * ——- PLANE ROTATION —— * SROT * 
x ® 2 & 
REEREEVTEEN REKEEKKEEE 
PURPOSE: To perform two dimensional rotations. 


CALL FORMAT: CALL SROT(N,SX,I,SY,J,C,S) 


PARAMETERS : N = Integer count of elements in SX and SY 
SX = Floating-point input vector of first components 
= (On output) first components of rotated vector 
I = Integer step increment for SX 
SY = Floating-point input vector of second components 
= (On output) second components of rotated vector 
= Integer step increment for SY 
Floating-point input scalar cosine 
= Floating-point input scalar sine 


nna 
iT} 


DESCRIPTION: SX(m) = C*SX(m)+S*SY(m) 
SY¥(m) =-S*SX(m)+C*SY(m); for m=l1 to N 


EXAMPLE: 


f 
Ne 


vs 
SX( INPUT) 
SY( INPUT) 
SX(OUTPUT) : 
SY(OUTPUT) : - 


ee ee oe 
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BRAEEEERARER RERRALEEEAE 
& * & 
* SROTM * -—- MODIFIED GIVENS ROTATIONS ——— * SROTM * 
* 2 & 2 
RRRARERAKERE RRERKREKREE 
PURPOSE: To perform two-dimensional rotations using 


the rotation matrix constructed froma 
parameter vector according to the modified 
Givens scheme. 


CALL FORMAT: CALL SROTM(N,SX,INCX,SY,INCY, PARAM) 


PARAMETERS : N = Integer element count 
SX = Floating-point input/output vector 
of first components 
INCX = Integer element step for SX 
SY = Floating-point input/output vector 
of second components 
INCY = Integer element step for SY 
PARAM = Five element floating-point input vector 
used to construct the rotation matrix 
H = H1l H12 
H21 H22. 


DESCRIPTION: SX(m) H11*SX(m) + H12*SY(m) 
SY¥(m) = H21*SX(m) + H22*SY(m), for m=l to N, where 


All, H12, H21, H22 = 
PARAM(2), 1.2, -1.9, PARAM(5) Or 
1.9, PARAM(4), PARAM(3), 1. or 


PARAM(2), PARAM(4), PARAM(3), PARAM(5) according to 
whether PARAM(1) = 1.9 or 9.9 or -1.M, respectively. 


If PARAM(1) is not equal to zero, one, or minus one, 
the routine returns with no action performed. This 
is equivalent to having the identity matrix as the 
rotation matrix. , 


EXAMPLE: 
N= 5 
SX(input) : @.8 1.28 -2.8 2.8 -4.9 
S¥(input) : 98.8 98.8 2.8 -2.8 -2.9 
PARAM : “1.8 1.89 71.8 1.8 1.8 
SX(output) : 9.9 1.09 8.80 8.8 -6.9 
SY¥(output) : 9.8 -1.0 4.8 -4.0 2.9 
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Rescaling continues until D1 and D2 are within the 
window. 

Output parameters PARAM(1,2,3,4,5) = 
(-1.9,H11,H21,H12,H22) and D1,D2,Bl are updated 
according to the scaling factors above. 


EXAMPLE: 
D1,D2,B1,B2 (input) : 4.988 3.9089 2.980 1.999 


D1,D2,B1 (output) : 3.368 2.526 2.375 
PARAM (output) : 9.8090 G.988 -8.580 9.375 G.9 
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REREEEKREARE RREKEARAREEE 
& « * * 
* SSWAP * ——— INTERCHANGES VECTORS —— * SSWAP * 
* & 2 z 
REERAKAKARIET REREREEKREX 
PURPOSE: To interchange elements of two real vectors. 


CALL FORMAT: CALL SSWAP(N,SX,1I,SY,J) 


PARAMETERS : N = Integer element count 
SX = Floating-point first vector for swap 
I = Integer element step for SX 
SY = Floating-point second vector for swap 
J = Integer element step for SY 


EXAMPLE: 
N = 3 
SX(INPUT) : 1.8 2.8 3.9 
SY(INPUT) : 9.8 8.8 7.9 
SX(OUTPUT) : 9.86 8.8 7.9 
SY(OUTPUT) : 1.86 2.8 3.98 
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RREEEEREKELR RREREEKRELR 


® . & x 
* ABP1 * --——- ADAMS-BASHFORTH PREDICTOR (ORDER 1) —- * ABP1l #* 
& x ®& * 
REEERERAER TPEETEEREER 
PURPOSE: To solve an initial value problem for a set of 


ordinary differential equations, using a first 
order predictor (Euler's) method. 
CALL FORMAT: _ CALL ABP1(N,H,Y,F,YP) 
PARAMETERS : N = Integer element count, number of equations 
H = Floating-point input scalar step size for t 
Y = Floating-point input vector of dependent 
variables Y(t) 
F = Floating-point i 
elements dY/dt =F 
YP = Floating-point output vector of predicted 
variables Y(t+H) 


DESCRIPTION: For the system of equations dY/dt=F(t,Y(t)), the 
SOlution at t'=t+tH is given by 


YP(m) = ¥(m) + H*F(m); for m=l to N 


This provides an explicit first order solution 

to the initial value problem for a given function 

at time t'=t+H, given the values of the function and 
its derivative at time t. The evaluation of the next 
derivative, corresponding to F(t+H,Y(t+H)) at the 
new time point, t'=t+2*H follows similarly. 


EXAMPLE: 
N = 3 
H = G.1 
x. : 18 2.8 3.8 
F : 1.0 1.9 1. 


YP" 5 1.1 2.1 3.1 
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® 2 


* ABP3 * 
* * 


RRREREREELR 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


' DESCRIPTION: 


EXAMPLE : 
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REREEREKRERE 
& * 


-——-— ADAMS-BASHFORTH PREDICTOR (ORDER 3) -——- * ABP3 #* 


® * 
RREERRREERE 


To solve an initial value problem for a set of 
ordinary differential equations, using Adams' third 
order predictor method. 


CALL ABP3(N,H,Y,F,F1,F2,YP) 


N = Integer element count, number of equations 

H = Floating-point input scalar step size for t 

Y = Floating-point input vector of dependent 
variables Y(t) 

F = Floating-point input vector of derivative 
elements dY/dt=F(t,Y(t)) 

Fl = Floating-point input vector of derivative 
functions at preceeding time tl=t-H 

F2 = Floating-point input vector of derivative 
functions at preceeding time t2=t-2H 

YP = Floating-point output vector of predicted 
variables Y¥(t+H) 


For the system of equations dY/dt=F(t,Y(t)), the 
solution at t'=t+H is given by 


YP(m) = Y(m) + (H/12)*(23*F(m)-16*F1(m)+5*F2(m) ); 
for m=l to N 


This provides an explicit third order solution 

to the initial value problem for a given function 

at time t'=t+H, given the values of the function and 
its derivative at t and its derivatives Fl and F2 at 
times tl=t-H and t2=t-2H, respectively. 

Evaluation of the next derivative, corresponding to 
F(t+H,Y(t+H)) at the new time point, t'=t+2*H 
follows similarly. 


N = 3 

H = 8.1 

Y 1.8 2.28 3.8 
iy 3 3.8 3.9 3.98 
Fl: 2.9 2.8 2.9 
F2 3: L.d 1.28 1.9 
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EXAMPLE: 


3.9 
3.0 
2.0 


3.9 
2.2 


3.9 


4.6 


4.9 


- 367 


Page A 
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DESCRIPTION: This routine integrates a set of N first order 
differential equations from t=A to t=B, given 
the initial values Y(t) and the values of the 
derivative functions dY/dt=F(t,Y¥(t)) calculated 
in the user supplied routine DFUNF(T,N,Y,F). The 
step size H is regulated to keep the maximum 
local error less than EPS. The maximum number of 
steps taken per call is limited by MAXIT. The 
maximum step size is limited by HMAX. Error 
return codes are provided to monitor the progress 
of the algorithm. 


REFERENCE: Burden,R.L., Faires,J.D., and Reynolds,A.C., 
"Numerical Analysis", Prindle, Weber & Schmidt, Inc., 
Boston, 1978: “Adams Variable Step-size Predictor- 
Corrector” Algorithm 6.5 


DFUNF (user supplied APFTN64 subroutine): 


SUBROUTINE DFUNF(T,N,Y,F) 


Cc 
Cc xxx DFUNF *** SAMPLE APFTN64 ROUTINE *** 
Cc 

DIMENSION Y(N), F(N) 
Cc 

DO 19 I=1,N 

F(I) = -Y¥(I) + T+ 1.9 
19 CONTINUE 

Cc 
Cc CORRESPONDS TO SOLUTIONS OF THE FORM 
Cc 
Cc Y(T) = YO * EXP(-T) + T 
Cc 

RETURN 

END 

INPUT: A = g.8 

B = 3.8 
N = 5 
HMAX = G.2 
MAXIT = 19¢ 
EPS = 1.9E-6 


Y(le1l), «er Y(5,1): 
1.8 2.8 3.8 4.89 5.8 
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RRAETREEAERZ . RREREEEREEE 
* eI 7 x 
* AMC1L * -—— ADAMS~-MOULTON CORRECTOR (ORDER 1) —— * AMC1l * 
& . 2 & *& 
RREREEAREE RKREKERERE 
PURPOSE: To solve an initial value problem for a set of 
ordinary differential equations, using a first order 


corrector (backward Euler) method. 


CALL FORMAT: CALL AMC1(N,H,Y,FP,YP) 


PARAMETERS ; N = Integer element count, number of equations 
H = Floating-point input scalar step size for t 
¥Y = Floating-point input vector of dependent 


variables Y(t) 

FP = Floating-point input vector of derivative 
elements dY/dt=F(t+H,Y(t+H) ) 

YP = Floating-point output vector of predicted 
variables Y(t+H) 


DESCRIPTION: For the system of equations dY/dt=F(t,Y¥(t)), the 
solution at t'=t+H is given by 


YP(m) = Y(m) + H*FP(m); for m=l1 to N 


This provides an implicit first order solution 

to the initial value problem for a given function 

at time t'=t+H, given the values of the function and 
its derivative at time t. The evaluation of the next 
derivative, corresponding to F(t+H,Y(t+H)) at the 


new time point, t'=t+2*H follows similarly. 
EXAMPLE: 

N = 3 

H = 9.1 

Y : 18 2.9 3.8 

FP: 18 1.8 1.2 


YP. :% 1.1 2.1 Jel 
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RRREKRARAEE RREKRERAARE 
4 2 & *& 
* AMC3 * ##—— ADAMS-MOULTON CORRECTOR (ORDER 3) --- * AMC3 * 
x & & x 
RREREEERER REAERRAEEEX 
PURPOSE: To solve an initial value problem for a set of 


ordinary differential equations, using Adams' third 
order corrector method. 


CALL FORMAT: CALL AMC3(N,H,Y,F,F1,FP,YP) 


PARAMETERS : N = Integer element count, number of equations 
H = Floating-point input scalar step size for t 
Y = Floating-point input vector of dependent 
variables Y(t) 

F = Ploating-point input vector o 
elements dY¥/dt=F(t,Y(t)) 

Fl = Floating-point input vector of derivative 
functions at preceeding time tl=t-H 

FP = Floating-point input vector of derivative 
functions estimated for t'=t+H 

YP = Floating-point output vector of predicted 
variables Y(t+H) 


+ 
4 


DESCRIPTION: For the system of equations dY/dt=F(t,Y(t)), the 
solution at t'=t+H is given by 


YP(m) = Y¥(m) + (H/12)*(8*F(m)-F1l(m)+5*FP(m)): 
for m=l1 to N 


This provides an implicit third order solution 

to the initial value problem for a given function 

at time t'=t+H, given the values of the function and 
its derivative at t, as well as, its derivatives at 
times tl=t-H and t'=t+H, corresponding to Fl and FP. 
Evaluation of the next derivative, corresponding to 
F(tt+tH,Y(t+H)) at the new time point, t'=t+2*H 
follows similarly. 


EXAMPLE: 

N = 3. 

H = 9.1 

Y : 1.9 2.4 3.8 
F : 2.0 2.08 2.8 
Fl se: 1.2 1.2 1.2 
FPP : 3.2 3.2 3.8 
YP os: Lee 2429 32525 
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EXAMPLE 


3.0 
3.8 
2.9 


1.9 
3.9 


3.9 
2.8 


4.9 


4.9 


= 375 


Page A 
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BRK(N,2) = 9.2 
and an input coordinate value x, BIN uses a binary 


l. The index IX that locates x within the 
coordinate value breakpoint table such that 


x(IX) <= x < x(IX+1) 


2. The product DR = D(IX) * R(IX) where 


i] 


D( IX) x(IX)-x 


1/(x( IX+1)-x(IX)) 


R(IX) 


When a program makes repeated calls to a breakpoint 
search routine (i.e., BIN or STEP), BIN should be used 
if it is suspected that the input coordinate x varies 
rapidly with respect to the values in the coordinate 
value breakpoint table. In this case, the binary 
(successive interval halving) search employed by 

BIN is more efficient than the step (nearest 

neighbor) search used by STEP. 


Refer to the function generation in Appendix E for 
additional information. 


EXAMPLE; 
N = 3 
BRK = 1.8 2.8 7.8 1.8 @.2 @.@ 
X = ek 
IX = 2 
DR = -9.92 


NOTE 
If x <= x(1) then IX = 1 
If x >= x(N) then IX = N-1l 
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DESCRIPTION: I(I+l), for I = 8 to N-l, is the value of the Ith 
modified Bessel functions of the first kind evaluated at 
the point X. Refer to equation 9.6.3 of Abramowitz 
and Stegun for the defining equation. 


K(I+l), for I = 9 to N-l, is the value of the Ith 
modified Bessel functions of the second kind evaluated at 
the point X. Refer to equation 9.6.4 of Abramowitz 

and Stegun for the defining equation. 


Warnings and errors are reported to the calling routine 
via IERR. If CBEIK completes normally, then IERR is 
set to zero. 


Warning condition codes are all between 1 and 99 
inclusive. The possible warning values and their 
meanings are as follows: 


PT atts ot 1 


DdD-=- WT Se hkAA = a aw meme Sw 
SVL GUNMVULaALLYII UL 


IERR = 1 N is too large 
outputs. In most instances, ABS(X) 
< 490.0; this means that the Nth order 
Outputs exceed the dynamic range of the 
machine. A suitable N is calculated, 
the Bessel function values are computed 
up to this new N, and the new N value 
is returned. 


Error condition codes are all greater than or equal to 
198. The possible error values and their meanings 
are as follows: 


IERR = 199 ISTEP and/or KSTEP are equal to 
-1,9, or l. 

IERR = 191 X does not lie within the boundary of 
(+/-699, +/-69@i). 

IERR = 192 N is equal to 1. N must be greater 
than or equal to 2. 


References: Abramowitz, M., and Stegun, I., “Handbook 
of Mathematical Functions", Ninth printing, 
pp. 358-369. 


Mason, J.P., “Cylindrical Bessel Functions 
for a Large Range of Complex Arguments", 
Computer Physics Communications, 39(1983), 
pp.1-ll. 
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REKREETEREE REEREKKERE 
r 4 & . x *= 
* CBEJYH * ———- COMPLEX BESSEL J, Y, AND H —— * CBEJYH * 
7 & s x 
RERARZAREE RERRERERKE 
PURPOSE: To compute the complex Bessel functions of integer 


order of the first kind, second kind, and one of 
the Hankel functions at a point X. 


CALL FORMAT: CALL CBEJYH (X, N, J, JSTEP, Y, YSTEP, H, HSTEP, IERR) 


PARAMETERS : X = Complex input scalar 

The point at which to evaluate all functions. 
This is restricted to the portion of the 
complex plane bounded by (+/-699,+/-6991i). It 
can take on the values (+/-6922, +/-69@i}. 

N = Integer input/output scalar 
On input, the number of function values to 
evaluate. If N <= 9g, then this routine returns 
with no action. If N = 1, then an error is 
reported. Note that the zero order function 
values are stored in the first elements of the 
complex output vectors. 
On output, the actual number of Bessel functions 
computed. The input value of N is modified only 
in the case where IERR = 1, if too many function 
values were requested. If IERR is not equal to 
1, then N is not modified on return to the 
calling routine. 

JSTEP = Integer input scalar 
Element step for J. This can be any value 
except -l, 9, or 1. This is the number of 
words to skip between complex elements. 

YSTEP = Integer input scalar 
Element step for Y. This can be any value 
except -l, 9, or 1. This is the number of 
words to skip between complex elements. 

HSTEP = Integer input scalar 
Element step for H. This can be any value 
except -l, @, or 1. This is the number of 
words to skip between compiex elements. 

J = Complex output vector 
The function values of functions @ through N-1l 
for Bessel functions of the first kind. 

Y = Complex output vector 
The function values of functions @ through N-1l 
for Bessel functions of the second kind. 

H = Complex output vector 
The function values of functions @ through N-1l 
for one of the Hankel functions. If the sign of 
the imaginary part of X is positive, then the 
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If the second Hankel function is desired when 
the imaginary part of X is nonnegative, it can be 
computed with the following equation: 


H2 = J-iY 


Similarly, the first Hankel function can be 
computed when the imaginary part of X is negative 


See ee oe ee ee 


Abramowitz, M., and Stegun, 


I., “Handbook 


of Mathematical Functions", Ninth printing, 


pp.358-368. 


Mason, J.P., 


for a 


"Cylindrical Bessel Functions 
Large Range of Compiex Arguments”, 


Computer Physics Communications, 39(1983), 


pp.1l-1ll. 


NNN Ww 


( 9.445474488934634E+999, 
(-8.657694535589279E+090, 
(-0.473368082G533007E+999, 


( 8.2274498948G4525E+9990, 
(-8.156496699689927E-Z01, 
(-8.535757979634719E+90G, 


g 


-9.496529947699122E+999), 
G.365928928827088E+99G), 
$.2473976415133086E+999) 


@.719158582GG15G5E+98D), 
G.6298G18G399G9G07E+9O9), 
8.577336957578681E+999 ) 


-8.519554586744886E-G91), 
-8.292666506762191E+G9S), 
-$.225978379G1979GE+GGG ) 
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EXAMPLE: 
See Appendix E for function generation. 
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F(x)=F(x(i))+(F(x(i+1) )-P(x(i)))*(x-x(1))/(x (itl) -x(1)) 


where 

x(i) = x-coordinate value at the i-th 
x-coordinate breakpoint 

x( itl) = x-coordinate value at the ({i+1)-th 
x-coordinate breakpoint 

x = Input x-coordinate value where the 
interpolated function value is desired 

F(x(i)) = Function value at x(i) 

F(x(itl)) = Function value at x(itl) 

F(x) = Interpolated function value at x 


and x(i) <= x < x(itl) 


EXAMPLE : 
See Appendix E for function generation. 
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desired functions, storing them in FVAL. Refer to the 
function generation in Appendix E for additional 


information. 


F(x)=F(x(i) )+(F(x(itl) )-F(x(i)))*(x-x(i))/(x(i+1)-x(i)) 


where 
x(i) 
x(i+l) 
F(x(i)) 


F(x(itl)) 
F(x) 


EXAMPLE: 


x-coordinate value at the i-th 
x~coordinate breakpoint 

x-coordinate value at the (i+l)-th 
x-coordinate breakpoint 

Input x-coordinate value where the 
interpolated function value is desired 


= Function value at x(i) 


Function value at x(i+1l) 


= Interpolated function value at x 


and x(i) <= x < x(i+l) 


See in Appendix E on function generation. 
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DESCRIPTION: FUN4 uses the indexes IX, IY, IZ and IW from the 
. breakpoint searches and the values NX, NY, NZ, and NW 

to find the first function value pairs in the function 
value breakpoint table. It then performs a linear 
interpolation between them by applying the formula 
given below eight times over the x-axis, four times 
over the y-axis, twice over the z-axis, and once 
over the w-axis. FUN4 repeats the process for all the 
desired functions, storing the computed function 
values in FVAL. Refer to the function 
generation in Appendix E for additional information. 


F(x)=F(x(i))+(F(x( itl) )-F(x(i)))*(x-x(1))/(x(it1)-x(i)) 


where 

x(i) = x-coordinate value at the i-th 
x-coordinate breakpoint 

x( itl} = x-coordinate value at the (it+tl}-th 
x-coordinate breakpoint 

x = Input x-coordinate value where the 
interpolated function value is desired 

F(x(i)) = Function value at x(i) 

F(x(itl)) = Function value at x(i+l) 

F(x) = Interpolated function value at x 


and x(i) <= x < x(i+tl) 


EXAMPLE: 
See Appendix E for function generation. 
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C.W.Gear, “Numerical Initial Value Problem in 
Ordinary Differential Equations", Prentice-Hall, 1971. 


RKGIL performs integration for given time, step 

size, and integration steps. The right-hand subroutine 
DFUN can be coded in either APFTN64 or APAL64. The 
parameter~passing method employed by RKGIL requires that 
DFUN be coded in APFTN64. As such, RKGIL relies on 
assumed procedure entry conventions, because APFTN64 
automatically generates code using this convention. 

If DFUN is written in APAL64, the user must resolve the 
parameters correctly. 


At output, vector V contains the numerical solutions 
while T8 contains the new value of the independent 
variable; i.e., TO=TH+M*H. 

poececne calls to RKGIL can cause eaewee proprens: 


and must take care speci eyida the H parameter. 


EXAMPLE : 
Solve the following second-order differential equation 
Y'* = ~-4.g*Y 


with initial conditions 


starting at TO = 9.9 with H = 9.1 for 32 iterations. 


An equivalent system of first-order differential equations 
can be written in the form 


DV(1) V(2) 


DV(2) -4.9*V(1) 


with initial conditions at the point 9.9% of 


V(1) 1. 


V(2) 0.0 
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REREKRERER RAREREAKEEE 
bY . & ® 
* RKGTF * ———- R-K-GILL-THOMPSON INTEG.(ORDER 4) —— * RKGTF * 
* * aot a * 
RRERERAREKE xkezxzkeaerzrzke 
PURPOSE: To solve an initial value problem for a set of 


ordinary differential equations, using the fourth 
order Runge-Kutta-Gill method as described by 
Thompson. 


CALL FORMAT: CALL RKGTF(T,N,Y,F,Q,H,M) 


PARAMETERS : T = Floating-point input scalar independent variable, 

initial value of t 

N = Integer input element count, number of equations, 
Gimension of Y, F and Q 

Y = Floating-point input/output vector of dependent 
variables (Y(t)) 

F = Floating-point working vector of derivative 
functions dY/dt=F(t,Y(t)) 

Q = Floating-point working vector used for 
temporary storage (must have length N) 

H = Floating-point input scalar step size for t 

M = Integer input scalar number of integration steps 
to be performed 


DESCRIPTION: For the system of equations dY/dt=F(t,Y(t)), the 


solution at each sten is given by 


aire Set Shit oS Sell VSS rn 


Y(m) = Y(m) 
+(H/6)*(k1+(2—sqrt(2))*k2+(2+sqrt(2))*k3+k4) 


for m=9 to N-1l, where 
k1l=F(T,Y) 


k2=F(T+H/2,Y+9.5*H*k1) 


k3 = F(T+H/2,Y+@.5*(-l+sqrt(2))*H*kl 
+9.5*( 2-sqrt(2))*H*k2) 
k4 = F(T+H,Y-@.5*sqrt(2)*H*k2 


+8.5*(2+sqrt(2))*H*k3) 


while the independent variable is advanced by H 
until T =T + M*H. 
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* ROT3 #* 
* * 


REREEEAAEE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


APPENDIX A 


REREEREEKRE 

* * 

--—- 3D ROTATION MATRIX, 3-ANGLE —— * ROT3 * 
x & 

RERREEELEE 


To form a three-dimentional rotation matrix as a 
product of three successive rotations about any three 
orthogonal axes. 


CALL ROT3(I,A,J,B,K,C,R) 


I = Integer input scalar axis indicator (plus or 
minus 1l=x, 2=y, 3=2) 

A = Floating-point input scalar angle(radians) of 
rotation about axis I 

J = Integer input scalar axis indicator (plus or 
minus 1=x, 2=y, 3=2) 

B = Floating-point input scalar angle(radians) of 
rotation about axis J 

K = Integer input scalar axis indicator (plus or 
minus l=x, 2=y, 3=2) 

C = Floating-point input scalar angle(radians) of 
rotation about axis K 

R = Floating-point output rotation matrix 
(3x3 matrix stored in column order) 


This routine calculates a 3x3 matrix as a product of 
three rotations about any three orthogonal axes: 


R(matrix) = R(K,C)xR(J,B)xR(I,A) 


1 i) g 


where R(1l,w) 
g cos(w) sin(w) 


8 -sin(w) cos(w) 


R(2,w) = cos(w) g -sin(w) 

g 1 g 

sin(w) cos(w) 

and R(3,w) = cos(w) sin(w) g 
-sin(w) cos(w) g 

8 g 1 
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REERERERER RREREERERK 
* & = & 
* SCSl1 * -——— SCALAR COS/SIN, TM INTERP.(ORD 1) -—-- * SCS1l * 
* a * *® 
RREREREREEE RREEEKEEEE 
PURPOSE: To rapidly calculate the cosine and sine of an 
angle(radians) using values stored in TMROM. 


CALL FORMAT: CALL SCS1(A,CA,SA) 


PARAMETERS: A = Floating-point input scalar angle(radians) 
CA = Floating-point output scalar cosine(A) 
SA = Floating-point output scalar sine(A) 


DESCRIPTION: CA = COS(A), SA = SIN(A) 


by interpolation of values stored in TMROM 

uSing a first order Taylor's series approximation. 
The returned values are accurate to approximately 
seven decimal digits. 


NOTE: For 15 decimal digits of accuracy at a slight 
decrease in speed, see the routine SINCOS. 


EXAMPLE: 
A = 1.2 
CA = 9.5493923 
SA = 9.8414719 
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An input coordinate value x, and the index IX from 
a previous call to STEP or BIN, STEP uses a step 
search to determine the following: 


l. The index IX that locates x within the 
coordinate value breakpoint table such that 


x(IX) <= x < x(IX+1) 


D(IX) = x(IX)-x 


R(IX) 


1/(x(IX+1)-x(IX)) 


When a program makes repeated calls to a breakpoint 
search routine (i.e., BIN or STEP), STEP should be 
used if it is suspected that the input coordinate 

x varies slowly with respect to the values in the 
coordinate value breakpoint table. STEP's 

nearest neighbor searching is more efficient than 
the binary (successive interval halving) search used 
by BIN. 


At tne outset, if no a priori knowledge of the value 
of x is available, the first call to STEP should 

set IX = N/2. An alternative strategy is to 

make the first call to BIN, which initializes 

IX, and then make subsequent calls to STEP. 


fer t ion generation in Appendix E for 
additional information. 
EXAMPLE: 
N = 3 
BRK = 1.8 2.8 7.86 1.89 G8.2 G§.98 
X 21: 
IX = 2 
DR = -49.G2 


NOTE 
If x <= x(l) then IX = 1 


If x >= x(N) then IX 


ul 
a 

t 
| a 
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RERAREEARELE RERREEEEKE 

& & & & 

* CONNMO * ——— NMO WITH CONSTANT VELOCITY -——- * CONNMO * 

r ® x 7 

REREREEREAE RRREEREKRETE 

PURPOSE: To apply normal moveout (NMO), with constant 
velocity, to a seismic trace. 


CALL FORMAT: CALL CONNMO(D,N,X,V,SR,NNMO) 


PARAMETERS: D = Floating-point output vector of trace 
sample times. 
N = Integer input scalar; element count for D. 
X = Floating-point input scalar; offset distance 
in feet. 
Vv = Floating-point input scalar; velocity in feet. 
SR = Floating-point input scalar; sample rate (ms). 


: 


Integer output scalar; index of initial sample 
of zero-fill in destination trace. 


DESCRIPTION: The normal moveout computation is described 
in seismic signal processing references, 
such as: 


“Introduction to Geophysical Prospecting" 
Dobrin, M.B., 
McGraw-Hill, Inc 
New York, N.Y., 


pp. 291-254. 


r+ e 
Os 


"Geophysical Signal Analysis" 
Robinson, E.A and Treitel, S., 
Prentice-Hall, Inc., 

Englewood Cliffs, N.J., 1989, 
pp. 1-35. 


The square-root computation inherent in the 
process is accomplished with one iteration of 
the Newton-Raphson method. 


Using a normal moveout process as defined by X, 
V, and SR, destination trace D is filled with 
the times from which to interpolate the adjusted 
trace values. 


The initial sample value of zero-fill in the 
destination trace is returned in parameter NNMO. 


A value of N+l for NNMO indicates no zero-fill. 
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RREERERERE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS: 


DESCRIPTION: 
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RRRERARERER 

J x. 

——— RECURSIVE FILTER -—— * IIR3G * 
® & 

ReRARAZELREEE 

To perform a recursive digital filter with up to 39 

poles and 38 zeros. 

CALL IIR39@(A,I,B,C,K,N,NZ,NP) 

A = Floating-point input vector of length N+NZ. 
Contains the data to be filtered. It will be 
assumed that A is indexed from -NZ to N-l. 

I = Integer input scalar. 

Element step for vector A. 

B = Floating-point input vector of length NZ+NP+1. 
Contains the coefficients of the Filter. It will 
be assumed that B is indexed from 9 to NZ+NP. 
B(@) contains the scalar multiple coefficient, 
B(1l) to B(NZ) contain the coefficients of the 
zeros, and B(NZ+1) to B(NZ+NP) contain the 
coefficients of the poles. 

C = Floating-point input/output vector of length- 
N+NP. 
Contains the filtered data. It will be assumed 
that C is indexed from -NP to N-l. On input, 
C(-NP) to C(-1) contain the initial values. On 
output, the computed values are contained in 
C(%) to C(N-1). 

K = Integer input scalar. 
Element step for vector C. 

N = Integer input scalar. 
Element count. 

NZ = Integer input scalar. 
Number of zeros. 

NP = Integer input scalar. 
Number of poles. 

Performs a recursive (IIR - Infinite Impulse 

Response) digital filtering difference equation as 

follows: 

C(t) = Sum({ B(j ) * A(t-j), j = BG to NZ J 
- Sum{ B(m+NZ) * C(t-m), m = 1 to NP J 
for t = 8 to N-l 

where the dimensions of the arrays are A(-NZ:N-1}, 

B(@:NZ+NP), C(-NP:N-1). The second sum equals zero if 

NP = @. 
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RARREUTEEREE RRAUEAEAERKEEE 
& & 4 , 
* KSMLV * —— K-TH SMALLEST ELEMENT IN VECTOR --—— * KSMLV * 
® x 2 x 
RREREEEREE RERTEREEEE 
PURPOSE: To find the k-th smallest element of a vector. 


s 


CALL FORMAT: CALL KSMLV(A,N,K,W,C) 

PARAMETERS: A = Floating-point input vector 

N = Integer element count for A 

K = Order of the element to be selected; K=l 
will select the smallest element; K=N will 
select the largest element; K=INT((N+1)/2) 
will select the median element. 

W = Work-space vector; the size of the work 
Space must be equal to N 

C = Floating-point output scalar 


DESCRIPTION: C = k-th smallest element of A(m), ml to N. 


The k-th smallest element of the vector stored 
in Main Memory starting at location A is found 
using an application of the divide and conquer 
strategy. The algorithm implemented is as described 
by Aho, Hopcroft, and Ullman: THE DESIGN AND 

ANALYSIS OF COMPUTER ALGORITHMS, Addison-Wesley 


TATA 


f 
1974, Pp. Q7-909, Tha resultant alamant ie etn ad 


& 14% whe WS GLE Wee ale We S44 KA ad ad Ne ed 


into Main Memory at location C. The original 
contents of the input vector are lost. 


The speed of this routine is data dependent. 


EXAMPLE: 


> 


8 5.8 2.8 -1.09 3.9 -39.6 18.7 5.8 
g 
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NNMO 


18 2.8 3.8 4.8 5.8 6.8 7.8 8.8 9.8 19.9 


(input) 
3.9 6.8 9.9 12.8 15.9 18.9 21.8 24.9 27.8 39.9 


33.8 36.9 39.8 9.89 86.8 G.8 €H8 GG GH GG 


D: 


4.0 


7.8 5.5 


7.8 8.5 19.8 9.5 


Ds: (output) 
2.5 4.8 5.5 
25 


18 6.8 G.8 GO.8 O88 G.86 G8 O88 G.G 
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EXAMPLE: 


2% 
29 
14 


1.82 2.8 3.6 4.0 5.6 6.0 7.8 8.8 9.8 19.8 


8 6.8 9.89 12.8 15.8 18.8 21.0 24.8 27.9 348.9 


(input ) 


3 
33.8 36.9 39.0 9.09 GG GH GH O88 G8 G.8 


D 


4.8 


5.5 7.8 8.5 19.80 18.9 7.@ 5.5 
1.89 8.8 6.89 8.89 G86 G8 86.8 GH GB 


D: (output) 
2.5 4.9 
2.5 


- 411 
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RREETEKERTE RRRAEREKAEE 
* * * * 
* RESNMO * -——— RESIDUAL NORMAL MOVEOUT -—— * RESNMO * 
, & x 
REREARELKEEE SRAUEREKERER 
PURPOSE: To stretch or squeeze a seismic trace via 


linear interpolation. 


CALL FORMAT: CALL RESNMO(A, B, C, NI, SR, D, NO, NNMO) 


PARAMETERS : A = Floating-point input vector; source trace 
of sample values. 
B = Floating-point input vector of input 
control times (ms). 
Cc = Floating-point input vector of output 
control times (ms). 
NI = Integer element count for B and C. 
SR = Floating-point input scalar; sample rate (ms). 
D = Floating-point output trace vector 


of sample values. 

NO = Integer element count for D. 

NNMO = Integer output scalar; index of initial sample 
of zero-fill in destination trace D. 


DESCRIPTION: The normal moveout computation is described 
in seismic signal processing references, 
such as: 


"Introduction to Geophysical Prospecting" 
Dobrin, M.B., 

McGraw-Hill, Inc., 

New York, N.Y., 1976, 

pp. 291-254. 


"Geophysical Signal Analysis" 
Robinson, E.A and Treitel, S., 
Prentice-Hall, Inc., 

Englewood Cliffs, N.J., 1989, 
pp. 1-35. 


Using a stretching/squeezing function as defined 
by B, C, and SR, source trace C is converted into 
destination trace D. 


The initial sample value of zero-fill in the 
destination trace is returned in parameter NNMO. 
A value of N+l for NNMO indicates no zero-fill. 


The speed of this routine is data dependent. 
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REERKEAEERE RREEKKEERE 
. J 5 = ; J 
* TMCONV * ——- CONVOLUTION (CORRELATION) ——— * TMCONV * 
= bd 2 * 
REEKREEEREZR RERARERREKE 
PURPOSE: To perform a convolution or correlation 

operation on two vectors, with the operand 


in Main Memory and the operator in TM. 


CALL FORMAT: CALL TMCONV (A, ITMB, C, N, M) 


PARAMETERS : A = Floating-point input vector (operand) 
ITMB = Integer address of B in TM 
C = Floating-point output vector 
N = Integer element count for C 
M = Integer element count for B 


(Integer element count for A = N+M-1) 


DESCRIPTION: C(m) = SUM(A(mt+q-1)*B(q)); 
for g=l to M and m=1 to N. 


NOTE: For convolution, the elements of operator 
vector B must be stored in TM in reverse 
order. 


TMCONV performs either a correlation or a convolution 
operation between the (N+M-1)-element operand (trace) 


- - fle \ 
yector A and the M-element cperator (kernel) vector 8B. 


The N-element result vector is stored in C. The 
result vector C may overlay the operand A. Vectors A 
and C reside in main data; vector B is in TMRAM. 

B must be placed in TMRAM using MTMOV or another 
Table Memory Library routine before calling TMCONV. 


NOTE: TMCONV is superior to CONV for M greater than 
or equal to 128; otherwise, CONV is superior. 
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* vol * 
& ci 
RREARAEKAEE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 


EXAMPLE: 
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RAREEERAEEE 


x * 
—— VECTOR ZERO TRENDS -—— * val * 
: ® . 
RREREKEERE 
To produce an output vector of 9's and 1' 
based on zero trends in the input vector. 
CALL V#1(A,1I,B,J,N,NPTS) 
A = Floating-point input vector 
bi = Integer element step for A 
B -= Floating-point output vector 
J = Integer element step for B 
N = Integer element count for A and B 
NPTS = Number of points of source to be 
considered in creating a 
destination point 
B(m) = 9.8 if ( (A(m-NPTS+1) .EQ. 9.9) .AND. 
(A(m-NPTS+2) .EQ. 9.9) .AND. 
(A(m) .EQ. 8.8) ) 
B(m) = 1.8 otherwise . 
for m = NPTS to N 
(Note that B(1l) = ... = B(NPTS-1) = 1.9) 
The vector A is scanned. If the current point 
of A and the last NPTS-1 points of A are 9, then 
the current point of B is set to zero. Otherwise 


the current point of B is set to 1.9. The 
resultant vector B is useful in stacking operations. 


N 16 

NPTS = 3 

A: 18 2.8 G68 G98 5.8 8.8 GB G.B 
9.9 19.08 11.812.8 13.80 9.09 8.8 O.D 

B: 180 1.8 1.86 1.80 1.86 1.86 1.8 G.8 
9.860 1.8 18 1.86 1.86 1.8 1.80 6.8 
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A(1) should be equal to 9.8, and all other values 
of A(i) and B(i), for i = 1 to NC, should be greater ° 
than @.@. 


The initial sample value of zero-fill in the 
destination trace is returned in parameter NNMO. 
A value of N+l for NNMO indicates no zero-fill. 


Routine NMOLI (linear interpolation) or NMOQI 
(quadratic interpolation) is generally called 


pS —_~ ---+ SS ee SS ee ewe Se we Se ae ee ee ee 


subsequent to routine VARNMO. 


The speed of this routine is data dependent. 


EXAMPLE: 
Nc = 4 
N = 188 
SR= 3.6 
K = 198.9 


A: Q.8 75.86 189.8 299.9 
B: 5@80.8 6000.9 7THO8G.09 8590.0 


NNMO = 68 


D( 1) D( 2) D( 3) DC 4) DC 5) D( 6) DC 7) 
28.98 28.87 289.59 21.53 22.83 24.44 26.39 


D(65) D(66) D(67) D(68) D(69) D( 192) 
192.39 195.38 198.37 G.88 0.09 eae 8.98 
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* VSCAND * 
* * 
REERERERKE 


PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


EXAMPLE: 
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RRREKAREEE 
* & 
* VSCANZ * 


* x 


-—— VECTOR SCAN FOR ZEROS -——— 


RRREKERERREE 


To scan a source vector and record in 


runnina foral 


a Finatinn vartonr a yz y 
om H44ds Oey Nee Wt le ob ae 


a destination ¥wew ee — 


of the number of zeros encountered. 


CALL VSCAN9(A,B,N) 


A = Floating-point input vector 

B = Floating-point output vector 

N = Integer element count for A and B 
B(m) = number of @‘s in A(i) through A(m); 


for m= 1toN 

Scans the N values of the source vector A. 
Records the cumulative total of zero values 
in the N elements of vector B. The resultant 
vector B is useful as a mute finder. 


N = 29 

A: 18 1.8 8.0 8.8 18 G8 O.8 GH 1.8 1.8 
18 86.8 6.8 OH G89 GH 18 G8 F.8 1.8 

B: 89.80 9.8 1.86 2.86 2.80 3.8 4.0 5.8 5.8 5.8 
5.8 6.8 7.86 8.8 9.89 19.8 19.8 11.8 12.0 12.9 
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RREEEERRAE RUEERREETAE 
* * & 4 
* CSFR2 * -—-—- SPARSE COMPLEX SYMMETRIC FACTOR -——— * CSFR2 * 
x . ‘® *& 
RRREREREEE : RRRERERERE 
PURPOSE: To perform an LDL' factorization of a complex, 


symmetric matrix A, where A is sparse and is 
represented in packed form. 


CALL FORMAT: CALL CSFR2(N,NS,S,ICP,IRN, ZTOL,WRK, IERR) 


PARAMETERS : N = Integer input scalar 
Order of the matrix A (must be greater than 1) 
NS = Integer input scalar 


Number of sparse elements (i.e., nonzero and 
fill-in elements) in the lower triangle of A 
S = Complex input/output array of length NS 
On input, S contains the sparse elements of 
the lower triangle of A in column order. On 
output, S contains the superposition of L and 
D with the diagonal elements reciprocated. 
ICP = Integer input array of length N+l 
Contains pointers into S to the first sparse 
element of each column with ICP(N+l) = NS + 1 
IRN = Integer input array of length NS 
Contains the row numbers that correspond to 
the elements in S 
ZTOL = Floating-point input scalar 
Zero tolerance value 
WRK = Complex scratch vector of length N 
IERR = Integer output scalar 
Error code whose values are: 
@ - Normal termination 
1 - Routine aborted because a diagonal 
element was computed to be zero (i.e., 
its absolute value squared was less than 
or equal to ZTOL) 
2 ~ Routine aborted because N < 2 


DESCRIPTION: This routine factors A into LDL‘ where L 
is a lower triangular matrix with ones on its 
diagonal, D is a diagonal matrix, and L' is the 
transpose of L. The factorization is performed 
without any row or column interchanges. 
L and D are superpositioned by suppressing 
the ones on the diagonal of L; i.e., if the 
Superposition of L and D is denoted by C, then 
C=L+bD- 1. The sparse elements of the super- 
position of L and D are stored in the corresponding 
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Thus the superposition of L and D with the diagonal 


elements of D replaced by their reciprocals is 


(8.5,-%.5) 
(9.8, 8.) 
(2.8,-1.8) 
(9.0, 8.8) 
(9.8, 8.8) 
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(9.5, 8.5) 
(9.8, G.8) 
(1.8, 1.9) 
(9.8, 9.9) 


(0.2,-8.4) 


(9.8, 8.8) (-8.25,8.25) 


(0.8, 8.8) 


(8.8, 


1.2) 


(9.25,8.8) 
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DESCRIPTION: First CSFR2 is called to factor A into LDL‘ where L 
is a lower triangular matrix with ones on its 
diagonal, D is a diagonal matrix, and L' is the 
transpose of L. The factorization is performed 
without any row or column interchanges. 

L and D are superpositioned by suppressing 

the ones on the diagonal of L; i.e., if the 
superposition of L and D is denoted by C, then 
C=L+oD- 1. The sparse elements of the super- 
position of L and BD are stored in the corresponding 
locations of S with the diagonal elements of D 
replaced by their reciprocals. L and D may contain 
nonzero elements where A contains zero elements. 
Collectively called "fill-in", these zeros must be 
included in S as input sparse elements of A. Failure 
to properly provide for fill-in results in 
undetermined action by this routine. 


Next, CSSV2 is called to solve the system in three 
steps: 


(1) Solve Lz=b for z (forward elimination) 
(2) Solve Dy=z for y 
(3) Solve L'x=y for x (backward substitution) 


This routine supercedes CSFS and differs from it in 
two important respects. First, CSFS2 is much faster 
than CSFS. Second, CSFS2 does not check to ensure 
that fill-in has been provided for properly; whereas, 
CSFS does. 


The scratch parameter WRK is not used in the current 
release of this routine; however, it has been 
retained for compatibility with CSFS. Thus, a scalar 
may be used in place for a vector for WRK. 


For a more detailed discussion, refer to Appendix C. 


The execution time for this routine is data dependent. 


EXAMPLE: Let A be the complex, symmetric matrix 


(1.9, 1.9) (8.8, 8.8) (3.0, 1.9) (08.8, 0.8) (8.8, 8.8) 
(9.8, 9.9) (1.8,-1.9) (8.0, 6.8) (2.8, 8.8) (8.8, G.8) 
(3.8, 1.9) (8.8, 6.8) (8.8, 1.8) (08.9, 0.8) (8.9, O.f) 
(9.9, @.8) (2.0, 6.8) (8.0, 6.9) (08.9, 0.8) (2.9,-2.8) 
(9.0, 09.9) (8.8, ©€.0) (8.8, O.8) .(2.8,-2.8) (6.9, 2.8) 


OTE: It is known apriori that Fili-in occurs in 
element (4,4). 
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RERREERRER SREEARAEREEE 
& = & * 
* cssv2 * ——— SPARSE COMPLEX SYMMETRIC SOLVE -——— * cCSSV2 * 
= x > d * 
RREREERERE REE 
PURPOSE: To find the solution to the system Ax = b, where A is 


a sparse, complex, symmetric matrix that is LDL' 


factored and is represented in packed form. 


CALL FORMAT: CALL CSSV2(N,NS,S,ICP,IRN,BX) 


PARAMETERS : N = Integer input scalar 
Order of the matrix A (must be greater than 1) 
NS = Integer input scalar 


Number of sparse elements (i.e., nonzero and 
fill-in elements) in A 
S = Complex input array of length NS 
Contains the sparse elements of the super- 
position of L and D with the diagonal elements 
reciprocated. The elements are stored in 
column order. 
ICP = Integer input array of length N+1l 
Contains pointers into S to the first sparse 
element of each column with ICP(N+1) = NS +1 
IRN = Integer input array of length NS 
Contains the row numbers that correspond to 
the elements in S 


On input, BX contains the right-hand side 
vector b. On output, BX contains the solution 
vector x. 


DESCRIPTION: This routine solves the system Ax = b where A is a 
sparse, complex, symmetric matrix that is factored 
into LDL'. L is a lower triangular matrix with ones 
on its diagonal, D is a diagonal matrix, and L' is 
the transpose of L. L and D are superpositioned 
by suppressing the ones on the diagonal 
of L; i.e., if the superposition of L and D 
is denoted by C, thenC =L+D- TI. 


The solution process consists of three steps: 

(1) Solve Lz=b for z (forward elimination) 

(2) Solve Dy=z for y 

(3) Solve L'x=y for x (backward substitution) 
This routine supercedes CSSV. 
For a more detailed discussion, refer to Appendix C. 


The execution time for this routine is data dependent. 


FPS 869-—7482-991C Page A — 429 


APPENDIX A 


RREREREEAEE REKERERERK 
& > * & 
* CUFR2 * -——~ SPARSE COMPLEX UNSYMMETRIC FACTOR -——— * CUFR2 * 
*& & & x 
PkeL EL ELSE SF | RERERAEKREEE 
PURPOSE: To perform an LU factorization of a complex, 
unsymmetric matrix A, where A is sparse and is 


represented in packed form. 


CALL FORMAT: CALL CUFR2(N,NS,S,ICP,IRN, IDP,ZTOL,WRK,IERR) 


PARAMETERS : N = Integer input scalar 
Order of the matrix A (must be greater than 1) 
NS = Integer input scalar 


Number of sparse elements (i.e., nonzero and 
fill-in elements) in A 
S = Complex input/output array of length NS 
On input, S contains the sparse elements of A 
in column order. On output, S contains the 
Sparse elements of the superposition of L and 
U with the diagonal elements reciprocated. 
IcP = Integer input array of length Nel 
Contains pointers into S to the first sparse 
element of each column with ICP(N+1) = NS + 1 
IRN = Integer input array of length NS 
Contains the row numbers that correspond to 
the elements in § 
= Integer input array of 1 
Contains pointers into S to the diagonal 
elements 
ZTOL = Floating-point input scalar 
Zero tolerance value 
WRK = Complex scratch vector of length N 
IERR = Integer output scalar 
Error code whose values are: 
9 - Normal termination 
1 - Routine aborted because a diagonal 
element was computed to be zero (i.e., 
its absolute value squared was less than 
Or equal to ZTOL) 
2 — Routine aborted because N < 2 


hw ONT 


tJ 
hg 
) 
ct 
@ 


anc 
eng 


FPS 869-7482-S81C Page A - 431 


APPENDIX A 


The output parameters are: 


Ss = 9.5, -0.5, 3.8, 1.8, 8.5, 0.5, 2.8, 8.0, 
2.8, -1.8, 8.2, -G.4, 1.8, 1.8, -G.25, 9.25, 
2.8, -2.9, 6.9, 1.8, 9.25, 8.9 

IERR = @ 


Thus the superposition of L and U with the diagonal 
elements of L replaced by their reciprocals is 


(@.5,-8.5) (8.8, 0.8) (2.0,-1.9) (8.8, 8.8) (8.8, 8.8) 
(9.8, 9.8) (8.5, 8.5) (8.8, 8.8) (1.8, 1.8) (8.8, 0.2) 
(3.8, 1.8) (8.0, 0.8) (8.2,-8.4) (9.8, 6.8) (8.8, 8.2) 
(09.8, 0.8) (2.8, 0.8) (8.8, 8.8) (-G.25,8.25) (0.8, 1.8) 
(9.0, 0.8) (8.8, 0.8) (8.0, @.8) (2.8,-2.8) (8.25,0.8) 
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DESCRIPTION: First CUFR2 is called to factor A into LU where L is 
a lower triangular matrix and U is an upper 
triangular matrix with ones on its diagonal. The 
factorization is performed without any row or column 
interchanges. L and U are superpositioned 
by suppressing the ones on the diagonal of U; i.e., 
if the superposition of L and U is denoted by C, then 
C=L+U- 1. The sparse elements of the super- 
position of L and U are stored in the corresponding 


| bs fF CC extikh kh As = 1 
20CGCions Ci GS Witnh ene GiagGondas elements of L 


replaced by their reciprocals. L and U may contain 
nonzero elements where A contains zero elements. 
Collectively called "fill-in", these zeros must be 
included in S as input sparse elements of A. Failure 
to properly provide for fill-in results in 
undetermined action by this routine. 


Next, CUSV2 is called to solve the system in two 
steps: 


(1) Solve Ly=b for y (forward elimination) 
(2) Solve Ux=y for x (backward substitution) 


This routine supercedes CUFS and differs from it in 
two important respects. First, CUFS2 is much faster 
than CUFS. Second, CUFS2 does not check to ensure 
that fill-in has been provided for properly; whereas, 
CUFS does. 


For a more detailed discussion, refer to Appendix C. 


The execution time for this routine is data dependent. 


EXAMPLE: Let A be the complex matrix 


(1.8, 1.8) (8.0, 89.8) (3.8, 1.8) (8.8, 9.8) (8.8, 8.2) 
(9.8, O.8) (1.8,-1.8) (8.0, 9.8) (2.8, 9.8) (8.8, 8.8) 
(3.9, 1.8) (8.8, 8.8) (8.8, 1.8) (0.8, 08.8) (8.8, 8.8) 
(9.8, @.8) (2.0, O08) (8.8, 0.8) (08.8, 0.8) (2.8,-2.f) 
(@.0, O@.8) (G8.8, O.8) (8.8, O©.H8) (2.8,-2.8) (6.8, 2.8) 


NOTE: It is known apriori that fill-in occurs in 
element (4,4). 


Let b be the complex vector 
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PURPOSE: 


CALL FORMAT: 


PARAMETERS : 


DESCRIPTION: 
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REEKEREEKRE 
z x 
——- SPARSE COMPLEX UNSYMMETRIC SOLVE —-- * CUSV2 * 
* x 
RREERERERRE 


To find the solution to the system Ax 


caAmn 1 ex 
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factored and is represented in packed 


CALL CUSV2(N,NS,S,ICP,IRN,IDP,BX) 


N = Integer input scalar 
Order of the matrix A (must be greater than 1) 
NS = Integer input scalar 


Number of sparse elements (i.e., nonzero and 
fill-in elements) in A 
iS) = Complex input array of length NS 
Contains the sparse elements of the super- 
position of L and U with the diagonal elements 
reciprocated. The elements are stored in 
column order. 
Integer input array of length N+tl 
Contains pointers into S to the first sparse 
element of each column with ICP(N+1) = NS + 1 
Integer input array of length NS 
Contains the row numbers that correspond to 
the elements in § 
= Integer 1 
Contains pointers into S t 
elements 
BX = Complex input/output vector of length N 
On input, BX contains the right-hand side 
vector b. On output, BX contains the solution 
vector x. 


IcP = 


input asLay or 


This routine solves the system Ax = b where A is a 
sparse, complex matrix that is factored into LU. L 
is a lower triangular matrix and U is an upper 
triangular matrix with ones on its diagonal. L and U 
are superpositioned by suppressing the ones on the 
diagonal of U; i.e., if the superposition of L and U 
is denoted by C, thenC =L+U- TI. 


The solution process consists of two steps: 


This routine supercedes CUSV. 


For a more detailed discussion, refer to Appendix C. 
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J & ® x 
* RSFR2 * ——— SPARSE REAL SYMMETRIC FACTOR —— * RSFR2 * 
* & ; > 
REEKEEEEEEE REKTEEKREE 
PURPOSE: To perform an LDL' factorization of a real, symmetric 


Matrix A, where A is sparse and is represented in 
packed form. 


CALL FORMAT: CALL RSFR2(N,NS,S,ICP,IRN, ZTOL,WRK, IERR) 


PARAMETERS : N = Integer input scalar 
Order of the matrix A (must be greater than 1) 
NS = Integer input scalar 


Number of sparse elements (i.e., nonzero and 
fill-in elements) in the lower triangle of A 
S = Floating-point input/output array of length NS 
On input, S contains the sparse elements of 
the lower triangle of A in column order. On 
output, S contains the superposition of L and 
D with the diagonal elements reciprocated. 
ICP = Integer input array of length Ntl 
Contains pointers into S to the first sparse 
element of each column with ICP(N+1) = NS + 1 
IRN = Integer input array of length NS 
Contains the row numbers that correspond to 
the elements in §S 
ZTOL = Floating-point input scalar 
Zero tolerance value 
WRK = Floating-point scratch vector of length N 
IERR = Integer output scalar 
Error code whose values are: 
9 - Normal termination 
1 - Routine aborted because a diagonal 
element was computed to be zero (i.e., 
its absolute value was less than or 
equal to ZTOL) 
2 ~- Routine aborted because N < 2 
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Then the input parameters are: 


N = 18 
NS = 22 
s = 8.9, 8.9, 16.9, 16.08, 32.8, 88.8, 16.9, 


24.8, 16.8, 8.8, 24.9, 8.0, 4.0, 16.9, 
32.0, 16.8, 88.9, 48.0, 8.0, 4.08, 9.8, 
=1.29 


TCD) Sly “2, ty 6,).. 8; “12s 145. L7y 203: 225 23 
IRN = i, Z, 18, 3, 4, 4, 5, 5, 6, 8, 6, 

8, 19, 7, 8, 9, 8, 9, 18, 9, 18, 10 
ZTOL = 1.9E-6 


The output parameters are: 


s = 9.125, 8.125, 2.8, 6.8625, 2.8, 9.8625, 1.9, 
@.125, 2.8, 1.8, -@.125, 1.8, -@.5, 9.8625, 
2.8, 1.9, 9.9625, 8.5, 8.25, -G.9625, 8.125, 
-$.93125 

g 


IERR 


Thus the superposition of L and D with the diagonal 
elements of D replaced by their reciprocals is 


8.125 

@.8 6.125 

9.9 8.8 8.8625 

G.09 8.9 2.8 0.2625 

6.8 8.8 6.6 1.2 @.125 

G8 688  G§.6 8.8 2.8 -@.125 

G@.8 6.6 6.8 6.6 0.6 6.6 6.9625 

G8 G8 G.8 8.G 18 180 2.8 6.8625 

G08 8.8 8.8 G.G 9.8 8.8 1.9 8.5 -8.8625 

9.8 2.8 G.9 G.8 G@.8 -O.5 8.8 8.25 8.125 -@.83125 
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DESCRIPTION: First RSFR2 is called to factor A into LDL’ where L 
is a lower triangular matrix with ones on its 
diagonal, D is a diagonal matrix, and L' is the 
transpose of L. The factorization is performed 
without any row or column interchanges. 

L and D are superpositioned by suppressing 

the ones on the diagonal of L; i.e., if the 
superposition of L and D is denoted by C, then 
C=L+bD- 1. The sparse elements of the super- 
position of L and D are stored in the corresponding 
locations of S with the diagonal elements of D 
replaced by their reciprocals. L and D may contain 
nonzero elements where A contains zero elements. 
Collectively called "fill-in", these zeros must be 
included in S as input sparse elements of A. Failure 
to properly provide for fill-in results in 
undetermined action by this routine. 


Next, RSSV2 is called to solve the system in three 
steps: 


(1) Solve Lz=b for z (forward elimination) 
(2) Solve Dy=z for y 
(3) Solve L'x=y for x (backward substitution) 


This routine supercedes RSFS and differs from it in 
two important respects. First, RSFS2 is much faster 
than RSFS. Second, RSFS2 does not check to ensure 
that fill-in has been provided for properly; whereas, 
RSFS does. 


The scratch parameter WRK is not used in the current 
release of this routine; however, it has been 
retained for compatibility with RSFS. Thus, a scalar 
may be used in place for a vector for WRK. 


For a more detailed discussion, refer to Appendix C. 


The execution time for this routine is data dependent. 


EXAMPLE: Let A be the symmetric matrix 
8.8 8.89 8.8 O88 8.8 GH F.86 G86 FG GF. 
9.86 8.86 6.8 880 8.86 68 88 8.809 G.8 16.9 
0.8 98.8 16.8 32.0 6.89 O89 G8 8.8 8.8 6.8 
9.8 9.8 32.06 86.09 16.89 689 O80 GF.80 8.86 98.8 
9.9 98.8 8.8 16.9 24.9 16.89 GO.8 8.89 8.8 8.8 
8.89 9.86 O@.8 8.8 16.8 24.6 6.89 8.8 8.8 4.9 
2.8 9.8 2.8 2.82 2.82 2.8 16.2 32.2 16.82 2.2 
9.68 86.80 G8 G68 8.80 8.8 32.6 88.8 49.9 8.0 
6.86 8.8 G8 8.8 G.8 S.86 16.80 49.0 4.8 6.8 
9.8 16.08 G86 6.8 6.86 48 6.89 8.86 9G.8 -1.25 
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Thus the superposition of L and D with the diagonal 


elements of D replaced by their reciprocals is 


8.8625 
1.9 
g.0 


G.8 2.8 
G.8 
G.G 


G8 


9.125 
2.08 


-$.125 


8.8625 


8.5 
8.25 


2.8 
1.9 
9. 


1.9 
8.8 
“8.5 


8.8 6.0 0.8 1.9 


G.8 


-8.98625 


@.125 ~@.93125 


8.9 9.9 


O.2 


is 


and the solution vector, x, 
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EXAMPLE: Let A 


RN RNwWRRNRQRNRA’ 
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NMNMQ®E NM NNN N OO 
BDNMNAN’MMNNM © QQ 
NWN NM FY QM NY MN MN QN RK 


le 


Then 


be the symmetric matrix 

9.98 8.89 6.86 8.8 G86 GOGH F.8 G.8 
$9 8.09 $09 8.86 8.86 G8 G.8 16.8 
16.9 32.90 9.80 99.0 G86 G8 S08 99.8 
32.9 88.8 16.9 9.9 9.8 G.8 8.8 6.6 
9.9 16.9 24.9 16.9 8.8 8.6 8.8 98.8 
9.8 9.80 16.09 24.0 98.8 8.08 6.8 4.9 
G@.6 6.0 6.6 6.6 16.6 32.9 16.6 6.6 
G8 G.9 8.8 8.8 32.0 88.8 408.9 8.08 
9.89 88 8.6 O8.8 16.80 409.6 4.0 8.8 
9.8 6.08 G.8 4.60 9.6 8.86 @.8@ -1.25 


the superposition of L and D with the diagonal 


elements of D replaced by their reciprocals is 


N) 
mn 


ie) 
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N 
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9.8625 

1.g 9.125 

g.f8 2.9 -@.125 

G.8 g.8 G.8 9.8625 

0.8 1.9 1.8 2.0 8.8625 

8.2 G.9 Gg. 1.9 g.5 -9.9625 
2.0 G.8 -g.5 g.8 8.25 8.125 -%.93125 
be the vector 

the input parameters are: 

= 19 

= 22 


9.125, 8.125, 2.8, 8.8625, 2.0, 8.8625, 1.9, 

8.125, 2.0, 1.9, -9.125, 1.9, -@.5, 8.8625, 

2.8, 1.8, 8.9625, 9.5, 9.25, -@.8625, 9.125, 

~%.83125 

= i, 2Z, 4, 6, 8, ii, i4, i7, 28, 22, 23 

= 1, 2, 18, 3, 4, 4, 5, 5, 6, 8, 6, 
8, 19, 7, 8, 9, 8- 9, 18, 9, 18, 18 

= 24.9, 8.9, 96.8, 288.8, 288.8, 296.9, 112.8, 

392.8, 28.8, 52.9 


Page A —- 447 


APPENDIX A 


RRRERREEEE REREREEEESE 
& x x rd 
* RUFR2 * -——— SPARSE REAL UNSYMMETRIC FACTOR --- * RUFR2 * 
* * = & 
PeeeeE SESE & ReEkEEEEEKE 
PURPOSE: To perform an LU factorization of a real, unsymmetric 
matrix A, where A is sparse and is represented in 


packed form. 


CALL FORMAT: CALL RUFR2(N,NS,S,ICP,IRN, IDP,ZTOL,WRK,IERR) 


PARAMETERS : N = Integer input scalar 
Order of the matrix A (must be greater than 1) 
NS = Integer input scalar 


Number of sparse elements (i.e., nonzero and 
fili-~in elements) in A 

s = Floating-point input/output array of length NS 
On input, S contains the sparse elements of A 
in column order. On output, S contains the 
sparse elements of the superposition of L and 
U with the diagonal elements reciprocated. 

ICP = Integer input array of length Nel 
Contains pointers into S to the first sparse 
element of each column with ICP(N+l1) = NS + l 

IRN = Integer input array of length NS 

Contains the row numbers that correspond to 

the elements in §S 


Tntanar innit ar 
2s Ceyos tit pusd i=" 


Contains pointers 
elements 
ZTOL = Floating-point input scalar 
Zero tolerance value 
WRK = Floating-point scratch vector of length N 
IERR = Integer output scalar , 
Error code whose values are: 
8 - Normal termination 
1 - Routine aborted because a diagonal 
element was computed to be zero (i.e., 
its absolute value was less than or 
equal to ZTOL) 
2 - Routine aborted because N < 2 


attr 
i 


r- 
) 
hg 

in 


1 
zt as 
into S to the diagonal 
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Then the input parameters are: 
N = 19 
NS = 34 
S = 8.9, 8.8, 16.9, 16.9, 32.9, 32.9, 89.9, 
16.8, 16.9, 24.8, 16.8, 8.8, 16.9, 24.9, 
8.9, 4.8, 16.9, 32.9, 16.9, 8.9, 8.9, 
32.89, 889.9, 49.0, 8.9, 16.9, 49.9, 4.9, 
9.0, 16.8, 4.8, 8.8, 6.8, -1.25 
ICP = ly 2, dp 6; “9p 1S) L7e 28y. 265. 3G, 35 
IRN = 1, 2, 18, 3, 4, 3, 4, Sy 4, 5S, 6, 
8, 5, 6, 8, 18, 7, 84 QO, Sy 6, 7, 
8, 9, 18, 7, 8, 9, 18, 2, 6, 8, 9, 19 
IDP = 1, 2, 4, 7, 18, 14, 17, 23, 28, 34 
ZTOL = 1.9E-6 
The output parameters are: 
S = 9.125, 9.125, 16.9, 9.9625, 32.9, 2.8, 9.9625, 
16.9, 1.8, 8.125, 16.9, 8.9, 2.9, -@.125, 
-8.0, 4.8, 8.9625, 32.8, 16.8, 1.9, 1.9, 2.2, 
9.8625, 8.9, 4.9, 1.8, 0.5, -@.9625, -2.2, 
2.0, -8.5, 8.25, 98.125, -#.93125 
IERR = 9 


Thus the superposition of L and U with the diagonal 
elements of L replaced by their reciprocals is 
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~~ 


G.8 8.8 G09 G.8 g.9 8.8 8.8 
G.G 8.6 8.8 §.8 8.8 8.9 2.28 
2.8 8.6 6.6 9.9 G.8 8.8 G.B 
@.9625 1.09 8.8 G.9 G.8 8.8 0.9 
6.8 9.125 2.6 9.6 1.8 8.8 G.8 
9.0 16.8 -G.125 9.0 1.0 B.8 9.5 

g.0 @.6 8.8 98.8625 2.9 1.G 8.2 
g.8 8.8 -8.8 32.8 $.9625 8.5 G.25 
g.8 9.6 6.8 16.9 8.8 -8.8625 9.125 
8.8 8.8 4.6 8.2 4.0 -2.8 -@.83125 
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DESCRIPTION: First RUFR2 is called to factor A into LU where L is 
a lower triangular matrix and U is an upper 
triangular matrix with ones on its diagonal. The 
factorization is performed without any row or column 
interchanges. L and U are superpositioned 
by suppressing the ones on the diagonal of U; i.e., 
if the superposition of L and U is denoted by C, then 
Cc =L+U- I. The sparse elements of the super- 
position of L and U are stored in the corresponding 
locations of S with the diagonai elements of L 
replaced by their reciprocals. L and U may contain 
nonzero elements where A contains zero elements. 
Collectively called "fill-in", these zeros must be 
included in S as input sparse elements of A. Failure 
to properly provide for fill-in results in 
undetermined action by this routine. 


Next, RUSV2 is called to solve the system in two 
steps: 


(1) Solve Ly=b for y (forward elimination) 
(2) Solve Ux=y for x (backward substitution) 


This routine supercedes RUFS and differs from it in 
two important respects. First, RUFS2 is much faster 
than RUFS. Second, RUFS2 does not check to ensure 
that fill-in has been provided for properly; whereas, 
RUFS does. 


For a more detailed discussion, refer to Appendix C. 


The execution time for this routine is data dependent. 


EXAMPLE: Let A be the matrix 
8.2 9.8 8.2 G.2 gf g.0 0.0 G.0 8.8 9.0 
9.0 8.8 8.8 9.0 2.8 G.0 8.8 9. 0.9 16.9 
g. 9.89 16.8 32.2 G.8 8.0 G.0 8.0 g.0 g.8 
0 9.9 32.9 88.9 16.2 8.0 8.8 g.0 G.0 g.8 
8.8 9.9 9.9 16.89 24.9 16.9 .2 8.9 g. g.9 
9.0 g.0 9. 9.89 16.89 24.9 8.0 8.0 G.0 4.9 
8.2 8.8 0.0 9.2 0.8 @.9 16.8 32.9 16.9 8.8 
g.2 9.8 0.0 9.8 8.9 8.9 32.9 88.0 49.9 8.9 
g.0 G.2 9.0 8.0 8.0 9.9 16.8 49.9 4.0 9.0 
G.8 16.9 9.2 0.0 2.0 4.9 0.0 8.9 9.8 -1.25 


NOTE: It is known apriori that fill-in occurs in 
elements (19,9) and (9,19). 
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Thus the superposition of L and U with the diagonal 


elements of L replaced by their reciprocals is 


G8 


8.8 


G0 


9.125 §.8 9.9 
9.09 98.125 9.8 
G.8 8.8 

¥.8 


9.0 
G.0 


8.0 


9.8 6.8 8.8 
G.8 8. 
2.8 
G.125 2.9 
16.9 


8.9625 2.9 


G. g.G g.8 


8.8 


8.8 


9.9625 1.9 


16 a) 


9.8 32.6 


g.8 
-G.5 


G.G 


g.2 


mS 
1.2 


g.0 


-$.125 9.9 


8.8 
8.8 


0.8 
G.8 


9.2 


8.9625 9.5 


@.25 


~8.8 32.9 
G.8 16.8 


8.0 
G.G 
G8 


9.9 
G.8 
G.8 


9.6 8.8 
8.8 
8.8 


G.8 


~G.4625 8.125 


-2.9 


8.0 


8.8 
G.8 16.9 


G.G 


-.83125 


4.9 


8.8 


4.9 


is 


and the solution vector, x, 


SSPVHPygs 
NoAFrFAMHABAE!E W 


- 455 


Page A 


FPS 8698-7482-981C 


APPENDIX A 


The execution time for this routine is data dependent. 


Let A be the matrix 


EXAMPLE 


G8 
G@.8 16.9 


G.8 
G.0 
G.8 


9.8 
8.8 
8.0 
G.8 


8.8 
8.8 
8.8 
8.8 


2.8 
G.8 
8.8 
G.8 


8.8 
8.2 


8.8 


G.8 8.8 
g.8 
@.8 16.8 32.9 


8.08 


G.0 


G.8 
g.8 


8.8 


8.8 
8.9 


9.8 32.09 88.8 16.8 


4.9 
ZG 
8.8 


G.8 
=Leeo 


8.8 


8.0 


0.8 


0.8 16.8 24.9 


G.9 
9.8 
G.2 
9.8 


G.8 
9.0 
g.8 
G.8 


G.G 
g.G 
9.8 


0.8 


89.8 16.8 32.8 16.8 


G. 


32.80 88.8 49.9 


8.9 


8.8 
8.8 
8.9 


8.8 


4.8 
9.8 


89.8 16.8 48.0 


4.0 


2.9 


8.8 


8.28 


B.2 


G.G 


8.9 16.9 


Then the superposition of L and U with the diagonal 


elements of L replaced by their reciprocals is 


g.9 


2.8 
g.8 
9.9625 2.9 


0.2 
G@.8 32.9 


G.125 6.2 


2.08 
GB 


GG. 8.8 G.8 B.D 


G.8 


9.125 9.8 


G.G 


G.0 
9.8 
8.8 
8.8 
G.B 
G.G 


9.9 
8.8 


G.8 
9.125 2.8 


$.9625 1.9 


16.9 


0.8 
8.8 
=0=5 


1. 
1.9 


0.8 


-9.125 9.2 
G.8 


16.8 


G.8 


8.8 


1.9 


G.8625 2.8 


G.8 


8.8 
G.8 
g.G 
2.8 


8.25 


8.9625 9.5 


8.9 
4.9 


“8.0 32.9 
8.8 16. 


4.0 


8.9 
8.8 


0.8 
G.8 


-$.8625 9.125 


-2.0 


9. 
B.G 


0.8 


-%.93125 


9.8 


G.8 


@.8 16.9 


Ctor 


Let b be the ve 


24.9 


8.9 


96.9 
288.9 


289.9 


296.9 
112.9 


392.9 


28.8 


52 2 


Then the input parameters are: 


8.8625, 


8.125, 16.8, 9.9625, 32.8, 2.0, 


9.125, 


~8.2, 


-9@.8625, -2.9, 


1.9, 8.5, 
-9.5, 9.25, 9.125, 


4.9, 


G.8625, 8.8, 
2.8, 


-9.@3125 
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EREEKEKEEE RERERERERE 
& * * = 
* SDOTPR * ——— SPARSE DOT PRODUCT —— * SDOTPR * 
® x & x 
RERERERAAE RRRAKREEKEEE 
PURPOSE: To calculate the dot product of a column of A 


with another vector, B, given a reali, sparse matrix, A, 
that is in packed format. 


CALL FORMAT: CALL SDOTPR(M,NP1,NS,S,IRN,ICP,IC,B,J,C) 


PARAMETERS : M = Integer input scalar 
Number of rows in A. 
NP1 = Integer input scalar 

Number of columns in A plus one. 


NS = Integer input scalar 
Number of nonzero elements in A. 

s = Floating-point input array of length NS 
Contains the nonzero elements of A stored by 
columns. 


IRN = Integer input array of length NS 
Contains the row numbers (in A) that correspond 
to the nonzero elements in S. 

ICP = Integer input array of length NP1l 
Contains pointers to the elements in S that are 
the first nonzero elements in each column of A. 
ICP(NP1L) = NS + l. 


IC = Integer input scalar 
Number of the column in A that is to be used. 
B = Floating-point input vector of length M 
J = Integer input scalar 
Element step for B. 
Cc = Floating-point output scalar 


DESCRIPTION: C = Sum({ B(IRN(k)) * S(k); k=ICP(IC) to ICP(IC+1)-1 ] 


EXAMPLE: Let A 3: 
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RREKETEEEEE ; REREEEERER 
* * * * 
* SITSOL * ——— SPARSE ITERATIVE SOLVER -— * SITSOL * 
x & * & 
KRREEEKAEKELT RERERAKRARKE 
PURPOSE: To solve a real, sparse, linear system A * X = B, 


wnere A is in packed, row-order format. 


CALL FORMAT: CALL SITSOL(N,NS,S,ICN,IRP,B,W,ZTOL,NCUT,IFLG, 


X, ITER, IERR) 
PARAMETERS : N = Integer input scalar 
Order of A. 
NS = Integer input scalar 
Number of nonzero elements in A. 
S = Real input array of length NS 


Contains the nonzero elements of A stored in 
row order. 

ICN = Integer input array of length NS 
Contains the column numbers (in A) of the 
corresponding elements in S. 

IRP = Integer input array of length N+l 
Contains pointers to the first element of each 
row of A in §S with IRP(N+1) = NS+tl. 


B = Real input vector of length N 
Contains the right-hand side. 
W ="Real input scalar 


Over relaxation coefficient. If W= 1.8, then 
the Gauss-Seidel method is used to solve the 
system. Otherwise, the successive over 
relaxation (SOR) method is used with a 
coefficient of W. 
Real input scalar 
Zero tolerance value. The solution is 
considered to have converged when every 
element of X is within ZTOL of its value on 
the previous iteration. 
NCUT = Integer input scalar 
Iteration limit. The routine will return 
after NCUT iterations if the solution has not 
converged. 
Integer input scalar 
Input flag: 
@ - Normal input 
1 - X¥ contains an initial solution 
2 - The routine is being reentered to 
perform additional iterations and the 
vectors S, ICN, IRP, B, and X contain 
the values that they had on return from a 
previous call to SITSOL. 


ZTOL 


H 
. Wy 
to 
Q 
" 
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EXAMPLE: Given the linear system 


A >: 4.9 8.2 2.8 
GG 8.0 9.8 
3.8 G.8 8.9 
Z.8 ~-3.9 9.9 
G.0 g.2 5.9 
G.0 g.0 8.8 
8 0.8 G.G 
G8 0.8 g.9 

and 


B : 8.8 ~5.9 7.9 
then the inputs are 


N 
NS 


oe 
[oe] 


A * X = B, where 


4 


~ 


e 
e 


NRE He FH QR 
MMM MM MQ QM Y 


RQRRQNN ®© QW 
MRM MM NM 


se 
e 


MN RN EMM M 
MWR MM MM MQ 


e 
e 


NR ®W WN M MN PY 
NN MMM MM YQ 


18.9 31.8 -22.9 4.9 


3.8, 3.8, 
5.2, -2.8, 
1.9, 2.@, 
4, l, 
3, 4, 
8, 7, 


IRP ¢ l, 3, 5, 8, 12, 17, 


5 Py &.2, -5 Oo, 7.8, 
W = 1.9 

ZTOL = g.9891 

NCUT = 29 

IFLG = g 


and the outputs are 


», « bd 2.99289, -1.999¢, 
1.9922, 3.9999, 


8 
g 


: 
w 
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8.9088, 
2.9999, 


3; 5, 

5; 6, 

6, 8 
Zip. 22% 


1.9908, 
S. SIG 


BQ RRA os} 
QNWRM MM MN MN N ® 


e 
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Input: 


EXAMPLE 


NS 


5.8 6.89 8.0 4.8 8.9 


A: 


3.8 g.8 
9.6 6.09 G80 GH G8 


6.89 8.09 G.89 G9 G.8 


9.6 G.8 GO. 


Output 


NS 


4.0 3.9 


5.89 9.8 6.8 


S 


IN: 


IP: 
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EXAMPLE: Input: 
ITYPE = 1 
M = 4 
N = 3 
NS =5 


S: 5.8 6.8 3.9 2.8 4.9 


4 
2 
oo 

te 
bo 
Pas 
~~ 
ge 


A: 5.86.09 8.93.09 6.0 6.089 O08 G.G6 2.86 G.86 G.8 4.8 
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Then the input is 


W) 
Zz 


IRN 


5 


ICP 


Ic 


8.8 3.8 2.80 2.8 
3.8 -5.9 


-7.9 


“1.9 


2.2 


4.8 


NC 


Output 
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Then the input is 


M = 6 
NP1 = 8 
NS = 12 
s ele 2e ~docel. “Ae mae. Su “SZe- Se Lan o2y. 43. 


ICP’ s. 1 3 4 5 7 1g 11 13 


bl b2 b3 b4 b5 b6 
where bl to b6 are the existing values in B 


w 


Output: 


B : bl b2 5. b4 -2. 3. 
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Output: 
NS = 3 
IERR = g 


S - 1.5 1425-=4.375 


IEN 3; 2 7 1g 
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Please detach cards along perforations. 


Your comments will help us improve the quality and usefulness of our publications. Please fill 
out and return this form. (The mailing address is on the back.) 


Title of document: 

Your Name and Title: Date: 
Firm: Department: 

Address: 

City: State: Zip Code: 
Telephone Number: ( ) Extension: 


I used this manual... . I found this material. . . 


QJ as an introduction to the subject accurate 


Q) to instruct a class written clearly 
OY to learn operating procedures well illustrated 


coooce 
ooocaczZ 


QO as a reference manual well indexed 


a other 


Please indicate below, listing the pages, any errors you found in the manual. Also indicate if 
you would have liked more information about a certain subject. 


= 
‘aus 
© 
es 
o 
ae 
bated 
= Q) as an aid for advanced training complete 
© 
O 
Y2 
ow 
fm] 
—_ 
<= 
fx] 
ae 


ARRAY is an independent society of people who use FPS products. Membership is free 
and includes a quarterly newsletter. There is an annual conference, as well as other 
activities. If you are interested in becoming an ARRAY member, please fill out and 
return this form. (The mailing address is on the back.) 


Your Name and Title: Date: 
Firm: Department: 
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City: State: Zip Code: 
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